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Abstract  of  the  Dissertation 

Combustion  and  Magnetohydrodynamic  Processes 
in  Advanced  Pulse  Detonation  Rocket  Engines 


by 


Lord  Kahil  Cole 

Doctor  of  Philosophy  in  Aerospace  Engineering 
University  of  California,  Los  Angeles,  2012 
Professor  Ann  Karagozian,  Chair 


A  number  of  promising  alternative  rocket  propulsion  concepts  have  been  devel¬ 
oped  over  the  past  two  decades  that  take  advantage  of  unsteady  combustion  waves 
in  order  to  produce  thrust.  These  concepts  include  the  Pulse  Detonation  Rocket  En¬ 
gine  (PDRE),  in  which  repetitive  ignition,  propagation,  and  reflection  of  detonations 
and  shocks  can  create  a  high  pressure  chamber  from  which  gases  may  be  exhausted 
in  a  controlled  manner.  The  Pulse  Detonation  Rocket  Induced  Magnetohydrody¬ 
namic  Ejector  (PDRIME)  is  a  modification  of  the  basic  PDRE  concept,  developed 
by  Cambier  (1998),  which  has  the  potential  for  performance  improvements  based 
on  magnetohydrodynamic  (MHD)  thrust  augmentation.  The  PDRIME  has  the  ad¬ 
vantage  of  both  low  combustion  chamber  seeding  pressure,  per  the  PDRE  concept, 
and  efficient  energy  distribution  in  the  system,  per  the  rocket-induced  MHD  ejector 
(RIME)  concept  of  Cole,  et  al.  (1995). 

In  the  initial  part  of  this  thesis,  we  explore  flow  and  performance  characteristics  of 
different  configurations  of  the  PDRIME,  assuming  quasi-one-dimensional  transient 
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flow  and  global  representations  of  the  effects  of  MHD  phenomena  on  the  gas  dy¬ 
namics.  By  utilizing  high-order  accurate  solvers,  we  thus  are  able  to  investigate  the 
fundamental  physical  processes  associated  with  the  PDRIME  and  PDRE  concepts 
and  identify  potentially  promising  operating  regimes. 

In  the  second  part  of  this  investigation,  the  detailed  coupling  of  detonations  and 
electric  and  magnetic  fields  are  explored.  First,  a  one-dimensional  spark-ignited  det¬ 
onation  with  complex  reaction  kinetics  is  fully  evaluated  and  the  mechanisms  for 
the  different  instabilities  are  analyzed.  It  is  found  that  complex  kinetics  in  addition 
to  sufficient  spatial  resolution  are  required  to  be  able  to  quantify  high  frequency  as 
well  as  low  frequency  detonation  instability  modes.  Armed  with  this  quantitative 
understanding,  we  then  examine  the  interaction  of  a  propagating  detonation  and 
the  applied  MHD,  both  in  one-dimensional  and  two-dimensional  transient  simula¬ 
tions.  The  dynamics  of  the  detonation  are  found  to  be  affected  by  the  application  of 
magnetic  and  electric  fields.  We  find  that  the  regularity  of  one-dimensional  cesium- 
seeded  detonations  can  be  significantly  altered  by  reasonable  applied  magnetic  fields 
(Bz  <  8 T),  but  that  it  takes  a  stronger  applied  held  (Bz  >  16 T)  to  significantly  alter 
the  cellular  structure  and  detonation  velocity  of  a  two-dimensional  detonation  in  the 
time  in  which  these  phenomena  were  observed.  This  observation  is  likely  attributed 
to  the  additional  coupling  of  the  two-dimensional  detonation  with  the  transverse 
waves,  which  are  not  captured  in  the  one- dimensional  simulations.  Future  studies 
involving  full  ionization  kinetics  including  collisional-radiative  processes,  will  be  used 
to  examine  these  processes  in  further  detail. 
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CHAPTER  1 


Introduction 

1.1  Background  on  Detonation  Engine  Concepts 

The  chemical  rocket  can  be  considered  the  oldest  technical  development  in  jet  propul¬ 
sion.  In  a  solid  propellant  rocket,  for  example,  the  exit  plane  momentum  is  due  to 
the  flow  of  a  hot  gas  created  by  the  rapid  burning  of  solid  fuel  composed  of  a  mixture 
of  a  fuel  and  oxidizer.  Gun  powder  emerged  in  China  around  AD  850  as  the  result 
of  accidental  discovery  by  Chinese  alchemist.  For  centuries  after  this  discover,  little 
was  done  in  the  advancement  of  rocket  propulsion,  until  in  1903,  a  Russian  school 
teacher  by  the  name  of  Konstantin  Tsolikowskyfl]  published  the  paper  ‘The  Investi¬ 
gation  of  Outer  Space  by  Means  of  Reaction  Apparatus”.  In  this  paper,  Tsolikowsky 
postulated  that  man  could  escape  the  clutches  of  earth’s  gravity  with  rockets.  His 
calculations  led  him  to  the  idea  of  multi-staging.  Subsequently,  he  went  on  to  dis¬ 
cuss  the  use  of  liquid  oxygen  and  liquid  hydrogen  for  those  purposes.  Piggy-backing 
off  of  these  ideas,  the  American  physicist,  Robert  Goddard,  designed  constant  pres¬ 
sure  rocket  combustion  chamber  nozzles  and  propellant  feed  systems.  Goddard[2] 
went  on  to  lead  the  advancement  of  liquid  fueled  rockets  and  in  1919  published,  “A 
Method  of  Reaching  Extreme  Altitudes” ,  which  not  only  provided  the  mathematical 
analysis  for  achieving  high  altitudes,  but  also  to  reach  the  moon.  He  devoted  much 
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of  his  effort  to  thrust  chamber  development  and  to  the  turbo-machinery  needed  for 
pumping  the  liquid  propellants  [3] . 

Still  today,  almost  a  century  later,  the  modern  aerospace  community  finds  itself 
facing  the  same  challenges  as  the  founding  fathers  of  modern  rocketry.  These  chal¬ 
lenges  include,  but  are  not  limited  to,  developing  lighter  and  more  efficient  propulsion 
systems  and  more  efficient  multi-staging  techniques.  The  present  studies  focus  on 
exploring  potential  rocket  propulsion  systems  that  take  advantage  of  magnetohydro- 
dynamics(MHD)  phenomena. 

In  1998,  Dr.  Jean-Luc  Cambier  proposed  a  novel  combined  cycle  propulsive  con¬ 
cept,  the  Pulse  Detonation  Rocket-Induced  Magnetohydrodynamic  Ejector  (PDRIME)[4], 
The  PDRIME  is  one  of  many  MHD  thrust  augmentation  ideas  that  shows  promise 
for  application  in  advanced  propulsion  systems.  Taking  advantage  of  the  unsteady 
wave  engine  concept  of  the  constant  volume  Pulse  Detonation  Engine  (PDE),  the 
PDRIME  utilizes  temporal  periodic  energy  divergence  into  a  seeded  air  stream,  then 
MHD  acceleration  for  thrust  augmentation.  Because  of  the  nature  of  the  unsteady 
waves  in  the  PDE,  the  PDRIME  does  not  need  heavy  turbo- machinery  to  pump 
liquid  propellants.  With  the  elimination  of  heavy  turbo- machinery,  paired  with  the 
energy  augmentation  of  the  MHD  accelerator  in  the  bypass  air  stream,  the  PDRIME 
could  potentially  be  able  to  achieve  the  velocities  necessary  for  Single-Stage-to-Orbit 
(SSTO)  flight,  thus  breaking  with  the  expensive  tradition  of  complicated,  expensive 
multistage  propulsive  systems. 

Another  concept  developed  by  Cambier  in  [5]  is  the  ‘Magnetic  Piston’.  This 
concept  involves  energy  extraction  from  the  expansion  portion  of  the  nozzle,  energy 
reintroduction  into  the  combustion  chamber,  followed  by  acceleration  of  the  com- 
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bustion  products  from  the  combustion  chamber.  The  ‘Magnetic  Piston’  builds  upon 
the  PDE  and  the  PDRIME,  but  also  gains  some  advantages  of  traditional  constant 
pressure  rocket  propulsion  systems.  These  and  other  alternative  rocket  propulsion 
concepts  will  be  explored  in  this  dissertation.  The  present  chapter  provides  the 
technical  background  for  these  concepts. 

1.1.1  Engine  Impulse  and  Efficiency 

Before  we  can  compare  the  properties  and  advantages  of  various  propulsion  systems, 
we  must  first  have  a  quantitative  means  of  expressing  various  properties.  Impulse 
and  efficiency  are  the  properties  we  are  most  interested  in.  Impulse  is  defined  as  the 
integral  of  thrust,  J7,  over  a  given  time  period: 

1=1  P(r)dr  (1.1) 

Jo 

Efficiency  in  a  propulsion  system  is  quantified  by  the  specific  impulse,  Isp,  which  is 
defined  as  Isp  =  where  Mpg  is  the  weight  of  the  propellant  used.  A  typical 
rocket  engine  uses  a  converging-diverging  (Laval)  nozzle  to  convert  high  pressure 
and  temperature  propellant  into  thrust.  The  larger  the  area  ratio  (AR),  ratio  of 
nozzle  exit  area  to  nozzle  throat  area,  of  the  nozzle  the  faster  the  exit  velocity  of  the 
propellant,  ue,  and  lower  the  exit  pressure,  Pe.  The  thrust  generated  by  rockets  is 
typically  expressed  as: 

T  =  mue  +  (Pe  -  Patm)Ae  (1.2) 

where  m  is  the  mass  flux  of  gas  exiting  the  nozzle  and  Ae  is  the  area  of  the  nozzle 
exit  plane.  For  optimal  thrust,  Pe  is  equal  to  the  ambient  pressure,  Patm ,  when  this 
is  achieved,  a  nozzle  is  said  to  be  ‘perfectly’  expanded.  When  traveling  in  altitudes 

which  range  from  sea  level  to  the  edge  of  space  (the  Von  Karman  line,  or  100  km), 
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there  is  a  large  variation  in  the  ambient  pressures,  so  that  most  often  Pe  ^  Patm- 
These  variations  in  ambient  pressure  lead  to  significant  losses  in  nozzle  efficiency. 

1.1.2  Pulse  Detonation  Engines 

The  Pulse  Detonation  Engines  (PDE)  in  simple  terms  is  an  engine  that  utilizes  deto¬ 
nation  waves  to  derive  its  thrust.  The  PDE  operates  in  a  cycle,  shown  schematically 
in  Figure  1.1.  In  this  simple  configuration,  a  stoichiometric  mixture  of  reactants  is 
placed  in  a  long  tube  with  an  open  and  closed  end  (i.e. ,  thrust  wall).  The  mixture 
is  then  ignited  from  the  thrust  wall  (a).  This  results  in  a  shock  and  deflagration 
wave  (subsonic  flame)  quickly  coalescing  into  a  detonation  wave  (supersonic  flame  - 
see  section  1.4)  propagating  into  the  reactant  mixture  (b).  The  high  pressure  region 
behind  the  wave  imparts  force  on  the  thrust  wall.  When  the  detonation  wave  reaches 
the  open  end  of  the  tube  (c),  it  will  be  reflected  back  into  the  tube  as  an  expansion 
wave  (d).  The  expansion  wave  propagates  into  the  tube,  while  at  the  same  time 
purging  the  combustion  products  from  the  tube  (e).  The  expansion  wave  will  then 
reflect  off  of  the  thrust  wall,  and  at  that  time  the  lowered  pressure  will  draw  fresh 
reactants  into  the  tube  (f).  The  expansion  wave  will  propagate  into  the  tube,  then 
reflect  back  from  the  open  end  as  a  compression  wave  (h).  This  compression  wave 
will  propagate  through  the  re-introduced  reactants  (i),  after  which  it  will  reflect  of 
the  thrust  wall  as  a  shock  (j)  and  a  new  cycle  will  begin.  There  are  many  advantages 
to  the  PDE  over  conventional  rockets  propulsion  systems.  The  conventional  rocket 
engines  have  to  pump  reactants  at  very  high  pressures  into  a  combustion  chamber, 
requiring  heavy  turbo  machinery,  while  the  PDE  does  not  need  heavy  turbo  ma¬ 
chinery,  but  rather  naturally  introduces  the  reactants  into  the  combustion  chamber 
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at  much  lower  pressures.  The  simplicity  of  the  concept  is  quite  attractive,  and  the 
PDE’s  have  been  tested  extensively  over  the  years  [22],  PDE’s  have  even  been  tested 
as  the  sole  source  of  propulsion  on  an  experimental  aircraft,  the  Scaled  Composites 
Long  E-Z[24],  which  used  an  abundance  of  off-the-shelf  parts. 

1.1.3  AJAX  and  RIME  concepts  for  thrust  augmentation 

Ejectors  have  been  considered  for  years  as  a  viable  method  of  thrust  augmentation 
for  various  aerospace  propulsion  systems.  Ejectors  rely  on  the  transfer  of  energy  from 
one  stream  (primary)  to  another  stream  (secondary).  Higher  thrust  can  be  achieved 
if  the  primary  stream  has  a  high  specific  energy  and  the  secondary  has  a  high  mass 
flow  rate.  A  generic  Rocket  Induced  Magnetohydrodynamic  Ejector  (RIME)  [25] 
consist  of  3  parts,  each  consisting  of  variable  area  stream-tubes;  (1)  the  generator, 
(2)  the  accelerator,  and  (3)  the  mixer.  The  streams  can  be  described  by  different 
power  plants,  for  example  a  rocket  stream,  a  bypass  tube,  and  the  mixer.  In  the 
RIME  described  by  [25] ,  the  “rocket  stream”  serves  as  the  MHD  generator  while  the 
“mixer”  is  the  MHD  accelerator.  The  MHD  generator  transforms  the  internal  energy 
of  a  fluid  into  electrical  power.  The  fluid  itself  is  a  conductor  of  electricity,  the  motion 
of  this  fluid  through  a  magnetic  held  gives  rise  to  electromotive  force  (drag)  and  how 
current  in  accordance  with  Faraday’s  law  of  inductance.  An  MHD  accelerator  uses 
the  same  principles,  but  in  this  particular  case  electrical  energy  is  applied  to  the 
system  resulting  in  an  electromotive  force  (thrust).  In  the  AJAX  system,  energy  is 
diverted  from  the  inlet  how  via  MHD  generation,  it  is  then  re-applied  after  the  fluid 
passes  through  the  combustor  via  MHD  generation [26]. 
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1.2  Plasma  Flows 


1.2.1  The  Lorentz  Force 

The  various  advanced  rocket  engine  configurations  explored  and  reviewed  in  this  dis¬ 
sertation  utilize  electromotive  forces  (Lorentz  forces)  generated  by  a  moving  charge  to 
augment  thrust.  The  Lorentz  Force  is  defined  as  the  Coulombic  attractive/repulsive 
force  between  single  or  collection  of  charged  particle(s)  moving  through  an  electro¬ 
magnetic  field: 


(1.3) 


F  =  j  x  B 


where  F  is  the  Lorentz  force,  j  is  the  current  density,  and  B  is  the  applied  magnetic 
field.  The  current  density  is  defined  by  Ohm’s  law: 


j  =  cr(E  +  u  x  B) 


(1.4) 


where  a  is  the  conductivity,  E  is  the  electric  field,  and  u  is  the  fluid  velocity.  Ex¬ 
pressing  the  electric  field  using  Ohm’s  law  and  using  vector  identities,  the  total  rate 
of  energy  deposition  into  the  fluid  can  be  expressed  as  : 


(1.5) 


where  the  first  term  on  the  right  side  is  the  heating  of  the  fluid  (dissipation)  and 
the  second  term  on  the  right  side  is  the  mechanical  power  obtained  from  the  Lorentz 
force  (non-dissipative).  The  separation  of  the  different  forms  of  power  expended  in 


the  fluid  becomes  important  as  we  evaluate  efficiency  later  in  this  document  [4], 
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1.2.2  Plasmas  and  Cesium  Ionization 

A  plasma  is  defined  as  a  collection  of  charged  particles  where  the  long-range  Conlom- 
bic  force  is  a  factor  in  determining  the  statistical  properties,  but  where  the  collection 
of  particles  is  low  enough  in  density  so  that  the  forces  exerted  by  a  particles  nearest- 
neighbor  is  less  than  the  long-range  Coulombic  force  exerted  by  the  particle’s  many 
neighbors.  Thus,  the  study  of  plasma  often  coincides  with  the  study  of  low-density 
ionized  gases.  An  ion  is  a  molecule  or  atom  which  acquires  enough  energy  to  liberate 
a  valence,  outer  shell,  electron  from  the  respective  molecule/atom.  A  collection  of 
ions,  liberated  electrons,  and  neutral  particle’s  form  a  plasma.  In  a  weakly  ionized 
plasma,  which  we  shall  investigate  throughout  the  remainder  of  this  study,  the  ion 
remains  in  close  proximity  with  a  liberated  electron,  thus  the  plasma  as  a  whole  can 
be  thought  of  as  charge  neutral,  but  locally  charged  [27]. 

One  of  the  necessary  properties  of  the  configurations  discussed  earlier  in  this  doc¬ 
ument,  the  AJAX,  the  RIME,  and  other  advanced  configurations,  is  the  ability  for 
the  working  fluid  to  conduct  electricity.  In  order  for  a  fluid  to  conduct  electricity  it 
must  be  at  least  partially  ionized.  These  advanced  configurations  achieved  ioniza¬ 
tion  through  thermal  ionization.  Thermal  ionization  follows  mass  action  laws  like 
any  chemical  reaction.  The  heat  of  ionization,  when  expressed  in  Kelvin  is  referred 
to  as  the  characteristic  temperature  of  ionization,  0*  =  Hionizaiton/k,  where  k  is 
Boltzmann’s  constant.  Most  common  gases  and  combustion  products,  i.e. ,  air,  CO, 
C02,  and  noble  gases,  have  high  characteristic  temperatures,  so  they  do  not  ther¬ 
mally  ionize  until  temperatures  in  excess  of  4000K  are  reached.  However,  if  an  alkali 
metal,  which  has  a  low  characteristic  temperature  of  ionization,  is  added  in  small 
amounts  (on  the  order  of  1  part  in  100  or  less)  thermal  ionization  can  be  achieved 
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at  temperatures  as  low  as  2000K.  This  process,  referred  to  as  seeding,  changes  the 
working  fluid  into  a  plasma  and  allows  the  previously  mentioned  configurations  to 
conduct  electricity  under  realistic  operating  temperatures  [28]. 

A  plasma  will  often  take  on  different  characteristics,  depending  on  the  temper¬ 
ature,  conductivity,  and  mean  velocity,  to  name  a  few.  Before  we  delve  into  the 
various  regimes  in  which  the  plasma  being  studied  will  exist,  we  must  first  look  to 
Maxwell’s  equations  for  free  charge [29]: 


VxH  =f +j 

V  x  E  =-§ 

V-B  =0 


(1.6) 


^  _  g  _  e(nj—ne) 

£0 

where  ne  is  the  number  density  of  electrons,  eo  is  the  permittivity  of  free  space,  H 
is  the  magnetic  field,  D  is  the  electric  displacement  held,  and  e  is  the  charge  of  an 
electron.  As  an  example,  if  a  one- dimensional  gas  were  to  how  perpendicularly  to  a 
magnetic  held,  the  induced  electric  held  could  be  expressed  as  Ey  «  \ uxBz ,  which 
can  be  shown  as  follows  for  the  idealize  case,  j/cr  — >  0: 


j  =  a  (E  +  u  x  B) 

E  «  — u  x  B  (1.7) 

Ey  UXBZ 


The  work  done  by  a  unit  volume  of  this  gas  moving  a  length,  L,  oriented  perpen¬ 
dicularly  to  the  magnetic  held,  B,  using  Equations  1.3  and  1.6,  would  produce  the 
following  [28]: 

FL  =  ^ aguB2L  (1.8) 
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The  energy  stored  in  a  magnetic  field  per  unit  volume  is  which  is  often 
referred  to  as  the  magnetic  pressure.  The  ratio  of  these  is 


Work  done  by  the  gas 
Energy  stored  in  the  field 


HvOgiiL  Rcm 


(1.9) 


The  Magnetic  Reynolds  number,  Rem,  is  a  good  measure  of  the  degree  in  which  a 
field  induced  by  gas  motion  compares  to  the  original  magnetic  field.  This  parameter 
plays  a  very  important  role  in  determining  the  performance  of  the  MHD  accelerator 
and  generator. 


1.3  Pulse  Detonation  Rocket  Induced  Magnetohydrodynamic 
Ejectors  and  other  Alternative  Configurations 

Spawned  from  the  concept  of  AJAX  [26]  and  RIME  [25],  Pulse  Detonation  Rocket 
Induced  Magnetohydrodynamic  Ejectors  (PDRIME)  and  other  alternative  configu¬ 
rations  are  meant  to  push  the  bounds  of  MHD  thrust  augmentation.  The  Pulse 
Detonation  Rocket  Engine  (PDRE)  is  the  core  of  the  various  alternative  engine  con¬ 
figuration  currently  under  review.  The  PDRE  is  composed  of  two  major  components: 
a  combustion  chamber  and  a  converging-diverging  nozzle  shown  in  Figure  1.2.  The 
combustion  chamber  introduces  reactants  on  the  front  end,  while  the  downstream 
end  connects  to  the  converging-diverging  nozzle.  The  converging  end  of  the  nozzle 
has  an  extremely  short  length  and  a  high  exit-to-throat  area  ratio,  AR  ~  16,  the 
significance  of  which  will  be  explained  momentarily.  A  typical  PDRE  cycle  is  at  the 
core  a  PDE  cycle,  where  reactants  are  introduced  to  the  combustion  chamber  and 
the  reactants  are  ignited  at  the  front  end.  A  detonation  wave  is  formed  and  propa¬ 
gates  downstream,  but  unlike  the  PDE  cycle  previously  mentioned,  the  detonation 
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wave  will  reflect  off  of  the  high  area  ratio  converging  section  nozzle  and  reflect  back 
as  a  shock  wave.  The  reflected  detonation  will  raise  the  chamber  temperature  and 
pressure  so  rapidly  that  the  expelled  gases  can  be  thought  of  as  started  from  this 
high  pressure  or  ‘blow  down’  state  [30].  The  fluid  is  then  accelerated  out  of  the 
nozzle  where  more  thrust  is  derived.  Over  a  cycle,  the  pressure  of  the  chamber  will 
decrease  as  more  products  are  being  expelled.  The  chamber  will  then  reach  a  critical 
pressure  in  which  reactants  will  be  introduced  at  relatively  low  pressure,  and  a  new 
cycle  will  commence. 

1.3.1  PDRIME 

The  PDRIME  is  actual  a  composite  of  some  of  the  systems  previously  discussed, 
that  is,  the  AJAX,  the  RIME,  and  PDRE.  The  PDRIME  is  physically  composed 
of  PDRE  (the  combustion  chamber  and  nozzle),  a  bypass  tube  that  sits  directly 
on  top  of  the  PDRE  and  magnets  which  are  placed  around  the  nozzle  as  well  as 
around  the  bypass  tube  which  is  illustrated  in  Figure  1.3  and  1.4.  The  effect  of  the 
strengths  of  these  magnetics  will  be  explored  later  in  this  dissertation.  In  a  PDRIME 
cycle,  reactants  as  well  as  a  gas  of  low  ionization  energy,  e.g.,  cesium,  are  introduced 
into  the  combustion  chamber  and  a  PDRE  cycle  will  commence.  The  fluid  in  the 
combustion  chamber  will  be  heated  sufficiently  to  ionize  the  seeded  cesium.  During 
this  process,  hot  products  and  ions  will  be  expelled  out  of  the  nozzle.  As  the  fluid 
expands  through  the  nozzle,  the  MHD  generator  in  the  nozzle  will  be  engaged.  The 
generator  will  extract  energy  from  the  fluid  moving  at  high  velocities  and  reduce  the 
Mach  number,  M  =  u/a ,  to  approximately  unity  at  the  nozzle  exit.  With  M  «  1 
at  the  nozzle  exit,  an  unsteady  shock  will  migrate  from  the  nozzle  exit  into  the 
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bypass  tube.  The  front  entrance  of  the  bypass  tube  is  seeded  with  cesium,  and  as 
the  unsteady  shock  migrates  upstream,  the  cesium  will  also  pass  through  the  shock 
and  be  ionized.  From  there,  the  MHD  accelerators  in  the  bypass  will  be  employed 
and  will  utilize  the  Lorentz  force  to  accelerate  the  fluid  out  of  the  bypass  tube.  The 
aim  of  this  configuration  is  to  take  energy  out  of  the  nozzle  and  more  effectively 
utilize  it  in  the  bypass  tube  in  order  to  increase  thrust  and  thus  gain  more  impulse 
and  efficiency. 

The  PDRE  has  many  good  performance  characteristics,  e.g.,  low  seeding  pressure, 
but  one  characteristic  that  we  wish  to  improve  is  the  nature  of  the  unsteady  pressure 
throughout  the  PDRE  cycle.  As  previously  described,  the  pressure  in  the  chamber  is 
quickly  increased  with  the  chemical  reactions  and  subsequent  reflection  of  detonation 
waves,  then  gas  is  expelled  from  the  combustion  chamber  and  expanded  through  the 
nozzle.  As  more  combustion  products  are  purged  from  the  combustion  chamber, 
the  pressure  in  the  combustion  chamber  drops  drastically  with  time.  In  a  PDRIME 
configuration,  this  drastic  drop  in  pressure  severely  handicaps  the  effectiveness  of 
the  MHD  accelerator  in  the  bypass  tube.  This  impairment  works  as  follows:  as  the 
pressure  drops  in  the  chamber  and  nozzle,  the  back  pressure  driving  the  nozzle  flow 
drops.  This  dropping  of  pressure  driving  the  nozzle  correlates  to  a  drop  in  pressure 
entering  the  bypass  tube  exit.  The  pressure  at  the  end  of  the  bypass  tube  supports 
the  unsteady  shock  wave  in  the  tube,  so  as  this  pressure  drops,  so  does  the  strength  of 
the  unsteady  shock  heating  the  fluid  in  the  bypass.  As  the  strength  of  the  unsteady 
shock  dies  down  is  strength,  the  temperature  jump  across  shock  is  reduced,  therefore 
less  ionization  takes  place.  Less  ionization  leads  to  lower  conductivity,  and  with 
lower  conductivity  the  MHD  accelerator  is  less  effective. 

In  order  to  prevent  the  negative  effects  of  the  unsteady  pressure  drop  through- 
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out  the  cycle,  MHD  forces  can  be  utilized  in  the  combustion  chamber  to  drive  the 
heated  fluid  out  of  the  combustion  chamber  at  a  more  constant  pressure.  This  con¬ 
figuration,  with  applications  of  energy  extracted  from  the  nozzle  into  the  chamber, 
for  a  “magnetic  chamber  piston”  configuration,  would  then  possess  the  advantages 
of  low  seeding  pressure  demonstrated  in  the  aforementioned  PDRE  [31]  as  well  as 
the  property  of  constant  pressure  possessed  by  a  conventional  rocket  engine  [5].  For 
the  basic  PDRIME  configuration,  this  combustion  chamber  has  chemical  reactions, 
ionization  of  seeding  ionizable  gas,  and  reflection  of  detonation  waves.  But  unlike  the 
PDRIME  cycle,  with  a  “chamber  piston”,  after  the  chamber  is  sufficiently  heated,  a 
MHD  accelerator  is  placed  around  the  combustion  chamber  to  create  the  piston  forces 
fluid  out  of  the  chamber  using  the  Lorentz  force.  The  ‘magnetic  piston’  concept  is 
illustrated  in  Figure  1.5,  and  the  PDRIME  with  the  ‘magnetic  piston’  is  illustrated 
in  Figure  1.6.  The  concept  of  a  ‘magnetic  piston’  was  first  introduced  by  Kolb  in 
1957[32],  In  his  experimentation  of  magnetic  shock  tubes,  Kolb  found  by  applying  a 
magnetic  field  to  a  plasma,  that  the  shock  waves  produced  were  stronger  than  that 
produced  in  the  absence  of  the  magnetic  field.  The  fundamental  interaction  of  a 
magnetic  field  with  shocks  and  detonations  will  be  explored  in  the  present  studies. 

1.3.2  Cambier’s  Quasi-ID  Model  and  Verification 

Cambier  has  performed  analysis  and  numerical  simulations  of  various  PDRE  configu¬ 
rations;  i.e. ,  the  PDRIME  and  the  PDRE  with  a  ‘magnetic  piston’. The  assumptions 
that  are  part  of  the  analysis  are  as  follows.  In  his  preliminary  studies,  Cambier 
selected  a  simple  configuration  where  the  fluid  velocity,  electric  field,  and  magnetic 
held  form  a  right  hand  coordinate  system  shown  in  Figure  1.7,  indicating  generator 
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(decelerator)  and  accelerator  configurations.  The  current  density  from  Equation  1.4 
can  now  be  expressed  as: 


jy  =  a(Ey  -  uxBz )  «  auxBz(Ky  -  1)  (1.10) 

where  K  is  the  loading  parameter,  i.e. ,  the  ratio  of  the  applied  electric  field  to  the 
induced  field:  Kv  =  E'!'  .  This  is  an  important  parameter  in  Cambier’s  simplified 
modeling  of  MHD  interaction  [4].  When  ux  >  0  and  0  <  Ky  <  1,  energy  is  ex¬ 
tracted  from  the  fluid,  the  Lorentz  force  is  negative,  as  is  jy,  and  the  device  acts 
as  a  “generator”,  shown  in  Figure  1.7(a).  When  Ky  >  1,  a  positive  application  of 
energy  takes  place,  the  Lorentz  force  is  positive,  as  is  jy,  and  the  device  acts  as  an 
“accelerator”,  shown  in  Figure  1.7(b).  It  is  the  generator  configuration  that  allows 
energy  extraction  as  shown  in  Figure  1.3,  with  energy  input  to  accelerate  the  flow, 
as  shown  in  Figures  1.4-1. 6.  In  the  present  studies  we  use  Ky  =  0.5  for  the  generator 
and  Ky  =  1.5  for  the  accelerator.  The  ideal  case  of  no  ohmic  heating,  j2/cx  «  0,  is  of 
particular  interest,  because  it  forms  a  very  simple  analytical  expression,  but  it  is  not 
always  valid.  In  order  to  neglect  the  ohmic  heating  as  compared  to  the  mechanical 
work,  the  following  condition  must  be  satisfied: 


j2/^ 

h 

Ey  ^ x 

=  \Ky-l\ 

u  •  (j  x  B) 

auxBz 

uxBz 

1  (1.11) 

Cambier  also  assumes  a  constant  magnetic  field  as  well  as  a  sufficiently  low  Rem. 
In  the  case  of  a  highly  conducting  plasma  (cr  1000  mhos/m),  which  can  be  seen  in 
the  combustion  chamber,  the  constant  magnetic  field  assumption  is  only  valid  if  the 
loading  factor  is  small,  which  also  leads  to  the  low  dissipation  approximation [4]. 
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1.4  Detonation  Waves  and  MHD  Effects 


The  foregoing  discussion  pertains  to  the  proposed  ability  of  magnetohydrodynamics 
(MHD)  to  affect  reactive  processes  via  seeding  the  flow  with  potentially  ionizing 
species  such  as  cesium.  While  this  notion  will  be  explored  in  a  global  sense  for  the 
PDRIME  and  alternative  detonation-based  engine  configurations  in  Chapter  5,  the 
ability  of  magnetic  and  electric  fields  to  influence  combustion  processes  requires  a 
more  detailed  examination,  in  particular,  including  the  effect  of  complex  reaction 
kinetics.  Hence  a  fundamental  understanding  of  transient  detonation  processes,  and 
the  ultimate  impact  of  MHD  on  these  processes,  is  required  and  will  also  be  examined 
in  this  thesis. 

The  study  of  detonation  waves  dates  back  to  the  late  19th  century,  where  Chapman [6] 
and  Jouguet[7]  modeled  detonations  as  a  shock  wave  supported  by  the  heat  release 
of  the  combustible  material  in  an  infinitely  thin  zone,  where  all  chemistry  and  dif¬ 
fusive  transport  takes  place.  Later  Zel’dovich[8],  Von  Neumann[9],  and  DoeringflO] 
independently  represented  the  detonation  as  a  confluence  of  a  one-dimensional  shock 
wave  moving  at  a  detonation  velocity,  followed  by  a  chemical  reaction  zone  of  finite 
length;  this  came  to  be  known  as  the  ZND  model  for  a  detonation  wave. 

While  the  true  structure  of  detonation  waves  inevitably  calls  for  representation  of 
multi- dimensional  effects  with  complex  reaction  kinetics,  the  simple  one-dimensional 
detonation  structure  provides  a  rich  spectrum  of  dynamical  features  which  are  worthy 
of  detailed  exploration  and  which  have  relevance  to  multi-dimensional  phenomena, 
e.g.,  cellular  detonations  [11].  Even  with  single  step  Arrhenius  kinetics  [12,  13,  14,  15], 
pulsations  or  instabilities  associated  with  a  ID  overdriven  ZND  detonation  may  be 
explored  in  detail,  with  important  physical  features  and  computational  requirements 
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established.  For  a  rapidly  initiated,  spark-induced  detonation,  where  the  detonation 
decays  from  an  over-driven  state  toward  the  self-sustaining  Chapman- Jouguet  limit 
cycle,  one  obtains  a  sequence  of  different  modes  of  physical  oscillation  between  the 
flame  and  shock  front.  The  numerical  analysis  of  this  effect  has  been  explored  previ¬ 
ously  by  Cambier[16]  using  complex  H2-&1T  detonation  kinetics  and  highly-resolved 
numerical  simulations,  but  with  only  a  spatially  and  temporally  second-order  con¬ 
vergence  rate  shock  capturing  scheme.  The  initiation  of  evolving  ID  detonation 
instability  modes  is  also  observed  in  calculations  by  Leung  et  al.  [17]  using  a  two-step 
chain-branching  reaction  model  [18]  and  a  Roe  scheme,  with  an  overdriven  ZND 
detonation  as  the  initial  condition.  Similar  calculations  with  a  second-order  accu¬ 
rate  slope-limited  centered  scheme  and  a  7-step  reduced  chemical  mechanism  for 
acetylene-oxygen  detonations  [19]  have  allowed  exploration  of  the  stabilizing  effect 
of  dilution  of  the  mixture  with  argon. 

In  the  present  study  (Chapters  6  and  7),  we  combine  higher-order  numerical 
methods  and  complex  reaction  kinetics  for  the  detailed  analysis  of  the  non-linear 
dynamics  associated  with  a  spark-induced  detonation.  While  simplified  one-step 
and  two-step  chemistry  models  have  provided  useful  guidance  in  elucidating  the  dy¬ 
namics  of  detonation  instability,  it  is  important  to  understand  the  influence  of  the 
complexity  of  a  realistic  reaction  mechanism,  since  energy  release  and  unsteadiness 
in  the  coupling  of  the  wave  front  and  induction  zone  can  affect  detonation  initiation 
or  failure  [20].  Moreover,  important  physical  processes  associated  with  deflagra¬ 
tion  to  detonation  transition  (DDT)  require  an  understanding  of  the  formation  and 
amplification  of  localized  explosion  centers  and  positive-feedback  flame  acceleration 
mechanisms  [21,  22,  23].  These  phenomena  are  often  easier  to  understand  through 
dimensionality  reduction,  while  preserving  some  of  the  complexity  of  the  physics  (i.e. 
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reaction  kinetics).  The  use  of  high-order  numerical  methods  can  also  be  a  powerful 
tool  in  the  analysis  of  such  complex  flows,  but  we  need  to  understand  the  interaction 
of  numerical  (spatial  accuracy)  and  physical  (chemical)  length  scales.  This  must  be 
done  before  adding  other  effects  such  as  species  diffusion  and  viscosity;  hence  our 
study  is  limited  to  reactive  Euler  flow.  Since  the  use  of  high-order  numerical  methods 
can  become  a  powerful  tool  in  studying  the  non-linear  detonation  dynamics,  it  is  also 
important  to  gain  a  good  understanding  of  the  effect  of  the  non-linear  algorithms 
on  the  flow  dynamics.  In  a  similar  fashion  to  many  previous  studies  of  detonation 
dynamics,  the  objective  of  the  present  work  is  not  to  provide  realistic  detonation 
simulations,  but  to  systematically  investigate  these  dynamics  through  the  addition 
of  increasing  complexity  in  the  models.  We  expect  of  course  that  adding  physical 
diffusion  will  eliminate  some  characteristic  length  scales  of  instability,  an  effect  which 
can  be  investigated  in  the  future. 

1.5  Goals  of  the  Present  Studies 

Among  the  goals  of  the  present  research  is  first  to  explore  alternative  configurations 
of  the  PDRIME  to  achieve  optimal  performance  from  the  MHD  augmentation.  The 
various  engine  configurations  used  in  our  research  took  on  a  multitude  of  forms,  these 
forms  are  shown  in  figures  1.2  through  1.6.  While  the  present  studies  described  in 
Chapter  5  focused  on  quasi-one  dimensional  simulations,  these  results,  together  with 
those  in  Zeineh  [33],  were  published  in  a  complete  study [34], 

Inherent  to  the  ability  of  the  PDRIME  and  its  modified  configurations  to  operate 
is  the  ability  of  an  applied  magnetic  held  to  affect  a  chemical  reaction.  Hence  in 
the  present  study  we  have  separately  studied  the  propagation  of  a  detonation  wave 
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with  complex  kinetics,  including  its  inherent  instabilities  (Chapter  6),  and  then  to 
examine  the  effects  of  the  cesium  seeding  and  an  applied  magnetic  field  on  the  dy¬ 
namics  of  the  detonation  (Chapter  7).  Thus  both  aspects  of  this  dissertation,  the 
simplified  modeling  and  the  detailed  detonation  simulations,  may  be  used  to  validate 
the  PDRIME  and  related  MHD  propulsion  concepts. 
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Figure  1.1:  The  generic  pulse  detonation  cycle,  (a)-(c)  represent  ignition  and  det¬ 
onation  wave  propagation,  (d)-(g)  represents  reflection  of  an  expansion  wave  from 
the  tube  opening  to  the  thrust  wall  and  back  to  the  tube  opening,  (h)-(j)  represents 
the  reflection  of  compression  waves  which  eventually  leads  to  the  re-ignition  of  the 
reactants  which  are  drawn  into  the  tube  at  stage  (f). 
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Figure  1.2:  PDRE  Pulse  Detonation  Rocket  Engine 


Figure  1.3:  Pulse  Detonation  Rocket  Engine  with  with  Nozzle  Generator  (NG)  with 
MHD  nozzle  generation  flight  configuration 
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Figure  1.4:  Pulse  Detonation  Rocket  Induced  MHD  Ejector  (PDRIME),  the  MHD 
accelerator  is  located  in  the  Bypass  Section 


Figure  1.5:  Pulse  Detonation  Rocket  Engine  with  Chamber  Magnetic  Piston  (CP) 
(from  Cambier[5]) 
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Figure  1.6:  Pulse  Detonation  Rocket  Induced  MHD  Ejector  with  Chamber  Magnetic 
Piston 
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Figure  1.7:  Lorentz  Force  (J  x  B)x  when  (a)  the  MHD  generator  is  on,  with  fluid 
moving  in  the  positive  x-direction,  applied  magnetic  field  in  the  positive  z-direction, 
and  current  flowing  in  the  negative  y-direction,  and  (b)  the  MHD  accelerator  is  on, 
with  fluid  moving  in  the  positive  x-direction,  applied  magnetic  field  in  the  positive 
z-direction,  and  current  flowing  in  the  positive  y-direction. 
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CHAPTER  2 


Governing  Equations  and  Physical  Phenomena 


An  accurate  description  of  the  governing  equations  is  presented  here,  describing  the 
flow  and  evolution  of  properties  associated  with  detonation  and  MHD  processes. 

2.1  Conservative  Formulation 

The  governing  equations  are  presented  here  in  the  differential  form,  but  later  will  be 
expressed  in  the  integral  form,  which  is  necessary  for  the  finite-volume  formulation 
used  in  some  simulations  of  the  PDRIME.  The  flow  equations  are  expressed  as  a 
hyperbolic  equation  with  a  source  term: 

^  +  Vn  •  F  =  O  (2.1) 

at 

where  the  Q,  F.  and  Cl  are  arrays  of  conserved  variables,  normal  component  of  the 
flux  density  of  Q,  and  source  terms,  respectively.  In  the  current  study  all  terms  on 
the  left  hand  side  (LHS)  of  Eqn  2.1  are  strictly  in  the  hyperbolic  form,  while  the 
right  hand  side  (RHS)  will  express  all  other  terms  which  will  be  referred  to  as  source 
terms.  Source  terms  can  describe  the  diffusion,  chemical  kinetics,  or  enforcement 
of  geometric  coordinate  constraint.  Operator  splitting,  which  will  be  discussed  in 

greater  detail  in  Section  3.1,  can  be  employed  to  compute  the  convective  contribution, 
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dt 


conv 


— V„  ■  F 


as  well  as  various  source  terms, 


(2.2) 


dQ 

dt 


mhd 


_  6 

“  Lmhdi  0 . 

at 


=  n 


kinetics 


kinetics 


(2.3) 


These  terms  then  can  be  combined  to  describe  the  total  change  of  the  conserved 
variables, 


9Q 

dt 


dQ 

dt 


dQ 

dQ 

dt 

conv 

+  dt 

mhd  Ui 

kinetics 


+ ... 


2.1.1  Single- Temperature(lT)  Hydrodynamic  Formulation 


(2.4) 


The  hydrodynamic  formulation  of  Eqn.  2.1  is  well  established  and  can  be  used  to 
describe  the  evolution  of  the  density,  velocity,  and  pressure  fields  of  the  fluid  where 
the  conserved  variables  and  normal  fluxes  can  be  described  as  such: 


Q  = 


/  p  \ 

rs 

pu 

v  E  / 


F  = 


\ 


Ps^n 

pu  ■  un  +  Pn 
y  (E  +  P)  un  J 


(2.5) 


where  the  total  mixture  density  p  =  n  is  an  arbitrary  direction,  un  =  u  •  n, 

and  the  total  energy  and  pressure  can  be  expressed  as: 


E  ^  ^  P s^int, s  T  „P U 


P=(7-l)(P--pn2 


(2.6) 

(2.7) 
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where  ejntjiS  is  the  specihc  internal  energy  of  the  sth  species  and  7  is  the  adiabatic 
index.  For  a  calorically  perfect  gas  (constant  cv)  the  specihc  internal  energy  can  be 
expressed  as: 

P  .  . 

& int,s  T  (2.8) 

7  - 1 

For  thermally  perfect  gases,  which  is  of  particular  interest  in  the  present  study,  the 
specihc  internal  energy  of  a  species  can  be  expressed  as: 

ef„M  =  J  cv,(T)dT  +  e0,,  (2.9) 

where  c„)S  is  the  specihc  heat  capacity  at  constant  pressure  for  the  sth  species  and 
the  eo,<i  is  the  specihc  internal  energy  of  formation  for  the  sth  species. 


2.1.2  Single- Temperature(lT)  Ideal  MHD  Formulation  ( Rem  — >  oo) 


By  combining  Maxwell’s  equations  and  the  induction  equation,  Equations  1.6  and  1.4 
respectively,  as  well  as  adding  a  zero  charge  separation  approximation,  ne  —  n*  ~  0, 
one  can  describe  the  time  evolution  of  the  magnetic  held  as: 


<9B 

~dt 


fi0a 


-VxVxB  — VxuxB 


convective 


diffusive 


(2.10) 


The  time  varying  B  contains  a  convective  term  which  behaves  as  a  hyperbolic  equa¬ 
tion  and  a  diffusive  term,  which  behaves  as  a  parabolic  equation.  Using  the  criteria 
set  in  previous  sections,  the  convective  term  will  be  treated  as  the  LHS,  of  the  form 
of  Eqn.  2.1,  and  the  diffusive  term  as  the  source  term.  If  the  conductivity  were  to 
be  extremely  large  such  that  o  — >•  oo,  then  the  diffusive  term  of  Equation  2.10  would 
become  zero,  leaving  only  the  convective  portion  of  the  equation.  When  this  partic¬ 
ular  case  is  cast  into  the  divergence  form  of  Eqn.  2.1  as  well  as  incorporating  the 
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magnetic  pressure  contribution  to  the  momentum  and  energy  equations,  i.e.  Lorentz 
force  and  Joule  heating  terms,  with  hydrodynamic  formulation  of  Equation  2.5,  it  is 
referred  to  as  the  ideal  MHD  formulation, 


/ 


P^n 


\ 


Q 


F 


puun  +  p* n  - 


n 


(2.11) 


B 


\E*  J  \  (E*  +  P*)un-±u-BBn 


where  Bn  —  B  n,  and  the  total  energy  and  pressure  can  be  expressed  as: 


Where  E*  is  the  total  energy,  P  is  the  mechanical  pressure,  and  the  total  pressure 
is  defined  as  P*  =  P  +  Pm.  Often,  is  referred  to  as  the  magnetic  pressure,  Pm. 
Without  the  presence  of  a  magnetic  field,  B  •  B  =  0,  the  formulation  of  Equation 
2.11  will  reduce  to  the  hydrodynamical  formulation  of  Equation  2.5. 

2.1.3  Two-Temperature(2T)  Formulation 

It  has  been  shown  by  Cambier[35]  that  under  certain  conditions  the  electrons  can 
be  heated  adiabatically  while  the  bulk  fluid  can  be  heated  non-isentropically.  Con¬ 
versely,  the  electrons  can  be  heated  non-adiabatically  while  the  bulk  fluid  is  heated 
adiabatically.  One  example  of  the  latter  is  when  a  microwave  is  used  to  excite  the 
electrons  in  the  fluid.  In  the  case  of  the  former,  a  fluid  may  pass  through  a  stationary 
shock  at  Mach  M  =  10,  but  while  the  electron  entropy  speed  is  the  same  as  that  of 


the  bulk  fluid,  the  speed  of  sound  of  the  electron  is  considerably  faster,  ce  = 
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where  the  ratio  of  molar  masses  of  Nitrogen  and  the  electron  is  ~  105.  In 

the  electron  reference  frame  there  is  therefore  no  shock,  Me  <C  1.  When  this  occurs, 
the  bulk  fluid  and  the  electrons  have  different  temperatures,  and  these  temperatures 
will  relax  on  time  scales  proportional  to  the  electron-heavy  particle  elastic  collision 
frequency.  The  two-temperature  MHD  formulation  (MHD2T),  Equation  2.1,  builds 
from  the  MHD  IT  formulation  but  contains  additional  terms  which  describe  the  evo¬ 
lution  of  the  electron  energy.  The  electron  thermal  energy  is  transported  as: 

dE 

+  V  (uEe)  =  — PeV  ■  u  (2.12) 


Because  the  electron  is  convected  at  u  and  |u|/ce  <C  0,  the  convection  of  the  electron 
is  subsonic  and  can  be  treated  isentropically.  This  allows  for  the  recasting  of  the 
electron  energy  into  the  electron  entropy,  Se,  which  is  a  conserved  quantity  and 
does  not  require  a  special  source  term  for  convective  transport,  where  the  conserved 
variables  and  normal  fluxes  can  be  described  for  the  two-temperature  formulation: 
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Q  = 


p  u 
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E* 


\  Se 


P^n 

puun  +  p*  n  -  Bn 
un  B  -  u  Bn 

(E*  +  P*)un-j-u-BBn 
k  Seun  J 


(2.13) 


where  the  electron  entropy  is  defined  Se  =  Es-  and  ye  =  |.  The  electron  energy,  total 
energy  and  total  pressure  are  defined: 


Seple  1  _  nekTe 

7e  -  1  7e  -  1 


v — -  1  9  B  •  B 

E  —  /  y  P spirit, s  +  2^U  ^ - 2 - ^  ^ 6 

s^e 
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Without  the  presences  of  a  magnetic  field,  B  B  =  0,  and  this  formulation  will  reduce 
to  the  two-temperature  hydrodynamical  formulation  (EULER2T). 


2.2  Overview  of  Source  Terms  for  Governing  Equations 


2.2.1  Combustion  and  Ionization  Reaction  Kinetics 


In  the  present  study,  the  kinetic  processes  which  include  the  chemical  reactions  of 
combustion  processes,  the  ionization  of  the  fluid,  and  the  temperature  relaxation  of 
the  electron  must  be  properly  captured.  In  the  case  of  chemical  reactions  and  ion¬ 
ization,  chemical  species  are  not  strictly  conserved,  but  particles  (chemical  elements) 
and  mass  are,  while  in  the  case  of  temperature  relaxation,  energy  can  be  transferred 
from  the  heavy  particles  to  the  electrons  and  vice-versa.  This  source  term  can  be 
represented  in  the  following  way, 


Q 


Ps 

cos 

pu 

^ kinetics 

0 

E 

lue 

Se  ) 

\  Se 

(2.16) 


where  us  is  the  production  of  the  sth  species,  uje  is  the  energy  production  due  to 
change  in  formation  energy,  and  Se  is  the  electron  entropy  production.  A  more 
detailed  description  of  these  terms  will  be  discussed  in  the  following  sections  as  well 
as  Section  2.3. 
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2.2.2  MHD  Transport 


2. 2. 2.1  Fixed  Magnetic  Field  ( Rm  ^0) 


In  the  MHD  source  formulation  that  is  used  to  study  the  PDRIME  configurations, 
and  some  of  the  MHD-detonation  interactions,  there  is  no  conservation  law  for  the 
magnetic  field,  but  rather  the  Lorentz  force  and  Joule  heating  are  incorporated  into  a 
hydrodynamic  formulation  of  the  governing  equations.  By  incorporating  the  Lorentz 
force  exerted  by  an  applied  magnetic  field,  Eqn.  1.3,  as  well  as  the  mechanical  power 
obtained  from  the  Lorentz  force  and  the  associated  Joule  heating,  Eqn.  1.5,  one  can 
recover  the  following  MHD  source  terms  for  a  fixed  magnetic  field: 


Q  = 


(  p  \ 

Ps 

pu 

\  E  ) 


( 


n 


mhd, fixed 


\ 


j  X  B 
J  E  ) 


(2.17) 


In  the  present  study,  when  this  approximation  of  MHD  is  used,  Cambier’s  [4]  sim¬ 
plified  MHD  model  will  be  implemented,  which  was  previously  discussed  in  Section 
1.3.2.  In  that  model,  the  system  includes  the  electric  and  magnetic  fields  as  orien¬ 
tated  in  Figure  1.7.  In  addition,  one  can  simplify  the  expression  for  current  density 
with  Equation  1.10,  such  that  the  x-component  of  the  Lorentz  force  becomes: 


(j  x  B)^  as  <mxB2JK,  -  1) 


(2.18) 


and  the  Joule  heating  term  becomes: 

j  •  E  «  au2xBl(Ky  -  1  )Ky  (2.19) 


This  formulation  will  be  employed  in  the  simplified  PDRIME  simulations  and  in 
some  of  the  detonation-MHD  studies. 
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2. 2. 2. 2  Resistive  MHD  (i?em  ~  0(1)) 


Seldom  do  real  problems  act  ideally,  as  described  in  the  perfectly  conducting  ideal 
MHD  approximation  of  Eqn.  2.11  or  the  perfectly  resistive  fixed  field  line  approx¬ 
imation  of  Equation  2.17.  Typically,  there  is  a  finite  conductivity  such  that  the 
diffusive  and  convective  terms  of  Eqn.  2.10  are  of  the  same  order  of  magnitude,  thus 
neither  can  be  neglected.  The  diffusion  of  magnetic  field  and  magnetic  energy  can 
be  expressed  as  the  following  source  term: 

^ mhd ,  diffuse 

where  the  evolution  of  the  magnetic  energy  is  prescribed  by 


V-  (  —  VB 

\/*0  <T 


(2.20) 


j-E  =  V-  (  — VT)  (2.21) 

\Po<?  J 

where  the  Maxwell  stress  tensor  is  defined  as  Ta/i  =  —  BO!BP .  A  more  detailed 

2/^0  mo 

discussion  of  magnetic  field  diffusion  will  be  presented  in  Section  2.4. 


2.3  Kinetics 

In  order  to  properly  resolve  the  chemical  reaction  and  ionization  processes,  one  must 
first  characterize  the  plasma  and  chemical  kinetics.  Using  detailed  balancing  one  can 
express  the  rate  of  species  production  and  destruction  in  the  following  manner. 

Us  =  ^2  Vrskfr  p]r°  -  ^2  "rshr  P^'3  (2.22) 

r  j  r  j 
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where  v”k  and  v'rk  are  the  coefficients  of  sth  species  in  the  rth  forward  and  backward 


reactions,  respectively,  and  kfr  and  kbr  and  the  forward  and  backward  chemical  rates 
of  the  rth  reaction.  When  the  temperature  of  the  fluid  is  near  thermal  equilibrium, 
detailed  balancing  can  be  used  to  determine  backward  rates: 


(2.23) 


where  A Gs  is  the  change  in  Gibbs  free  energy  of  the  sth  species.  In  cases  where  the 
temperature  is  far  from  equilibrium,  for  example  when  there  is  a  heavy  and  electron 
temperature,  detailed  balance  would  not  be  appropriate  to  determine  the  backward 
rates.  The  forward  reaction  rates  are  given  by  the  modified  Arrhenius  equation  of 
the  form: 


k  =  AT^expi-Q/T) 


(2.24) 


Where  A  is  the  Arrhenius  pre- factor,  rj  is  the  Arrhenius  coefficient,  and  0  is  the 
activation  temperature.  Using  the  form  of  Equation  2.16  with  Equation  2.22,  the 
energy  production  due  to  the  internal  energy  change  can  be  expressed  as  c he  = 

,s- 

2.3.1  Combustion  Kinetics 

2. 3. 1.1  Single  Step  Kinetics 

During  preliminary  testing,  we  simulated  a  single-step  H2  —  02  reaction,  reactants 
II)  and  02  form  product  H20 


H2  +  \o2  ->  H20 
35 
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The  H2O  production  rate  takes  the  Arrhenius  form  of  Equation  2.24  in  which  up  = 
Kpfe~e/T.  I11  this  particular  case,  pj  is  the  density  of  a  stoichiometric  mixture  of 
the  fuel,  H2,  and  oxidizer,  02,  and  pp  is  the  density  of  the  product,  H20.  The 
thermodynamic  properties  of  the  reactants  and  product  are  given  in  Table  2.1.  The 
differences  in  the  thermodynamic  properties  between  reactants  and  products  found 
in  this  table  are  due  to  the  differences  in  the  molar  mass,  M,  and  the  degrees  of 
freedom  associated  with  the  molecular  configuration  of  each  species. 

2. 3. 1.2  Detailed  Reaction  Kinetics 

In  the  present  study,  we  have  primarily  focused  on  detailed  kinetics  of  a  simple 
combustion  system  (H2  —  02).  Another  approach  to  the  single  step  reaction  in 
Equation  2.25,  commonly  chosen  in  fundamental  studies  of  detonation  dynamics,  is 
a  constant- volume  one-step  reaction  model,  in  which  the  entire  chemistry  is  described 
by  the  evolution  of  a  single  progress  variable  that  follows  an  exponential  relaxation 
with  a  characteristic  time-scale  given  by  the  induction  delay.  This  progress  variable 
is  also  associated  with  the  fractional  amount  of  heat  released  into  the  flow.  In  that 
model,  the  induction  delay  time,  rind,  follows  a  simple  exponential  fit,  tind  ~  e6a^T . 
The  delay  being  essentially  caused  by  the  need  for  a  sufficient  amount  of  radicals 
from  chain-branching  reactions,  and  the  production  of  those  being  an  endothermic 
process,  the  parameter  9a  in  this  formulation  is  an  averaged  activation  energy  of  the 
key  radical-producing  reactions. 

This  is  a  reasonable  approximation  to  the  chemistry  in  that  region,  albeit  within 
limits.  To  study  detonation  dynamics  more  completely,  we  have  used  the  detailed 
chemistry  to  compute  and  parametrize  the  induction  delay  as  a  function  of  initial 
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temperature  and  pressure,  with  the  mixture  held  hxed  to  stoichiometric  hydrogen- 
air.  As  in  a  previous  study  [16],  the  chemical  kinetics  of  a  dilute  hydrogen-air  mixture 
were  solved  here.  The  chemistry  includes  eight  reacting  species,  H2,  02,  H ,  O,  OH , 
H02 ,  H202 ,  H20,  and  the  non- reacting  diluent,  N2,  where  a  compilation  of  NASA 
and  JANNAF  thermo-chemical  data  is  obtained  from  [36].  As  prescribed  by  [37], 
38  elementary  reactions,  found  in  Appendix  A,  are  used  in  this  mechanism  and  the 
backward  rates  are  computed  from  equilibrium  constants.  As  shown  in  Figure  2.1,  the 
delay  does  follow  an  exponential  form,  f*  cx  a(P)e/3<'P^T ,  as  expected.  The  parameter 
/?  in  this  formulation  is  an  averaged  activation  energy  of  the  key  radical-producing 
reactions.  However,  this  approach  yields  unrealistic  profiles  of  the  post-shock  region, 
since  the  heat  release  is  gradual. 

A  better  description  is  obtained  with  a  two-step  reaction  model  [38,  39],  where 
the  heat  release  is  associated  with  a  second  progress  variable  whose  evolution  can 
start  only  at  the  end  of  the  induction  delay,  which  now  follows  a  linear  time  variation. 
While  this  two-step  model  allows  a  separation  between  the  induction  and  heat  release 
zones,  the  model  is  unsatisfactory  in  several  ways.  First,  the  rate  of  heat  release 
is  assumed  independent  of  temperature,  which  is  unrealistic,  as  the  flow  heating 
accelerates  the  combustion.  Generally  speaking,  a  stiff  differential  equation  for  the 
progress  variable  can  be  used  to  reproduce  this  non-liuear  effect,  but  the  dynamics 
can  be  different  from  the  real  conditions.  Second,  when  the  flame  is  accelerated 
towards  the  shock,  the  two  reaction  zones  (induction  and  flame)  start  to  merge, 
even  if  species  diffusion  is  neglected;  the  enforcing  of  two  separate  zones  with  a 
two-step  model  could  thus  modify  the  dynamics  of  the  strongly  coupled  shock-flame 
system.  It  is  for  these  reasons  we  will  utilize  detailed  reaction  kinetics  to  perform 
the  simulations  in  the  present  study,  per  Equations  2.22  -  2.24  and  the  full  H2  —  02 
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combustion  mechanism  described  in  Appendix  A. 

2.3.2  Ionization  Kinetics 

The  existence  of  plasma  in  a  working  fluid  is  studied  via  the  field  of  magnetohydro¬ 
dynamics.  If  a  working  fluid’s  characteristic  temperature  of  first  ionization,  0,;,  is 
high  relative  to  the  fluid  temperature,  it  is  imperative  to  seed  the  fluid  with  a  species 
with  a  sufficiently  low  characteristic  temperature  to  create  the  flow  of  ionized  gas. 
The  principle  components  of  the  working  fluids  in  this  study,  air,  H2l  02,  and  H20, 
have  high  characteristic  temperatures  relative  to  the  fluid  temperature  in  the  scope 
of  the  present  study.  A  sampling  of  ionization  temperatures  are  shown  in  Table 
2.2.  As  prescribed  by  Cambier  in  [40],  cesium  is  chosen  as  the  seeded  species  due 
to  its  low  characteristic  temperature  of  ionization.  In  the  modeling  of  the  ionization 
and  recombination  of  cesium  in  the  working  fluid,  we  start  with  a  simple  three-body 
reaction  mechanism: 

Cs  +  M  ^  Cs+  +  M  +  e~  (2.26) 

where  M  in  this  particular  case  represents  a  third-body  species.  In  the  present  study, 

it  is  extremely  important  to  be  able  to  calculate  the  conductivity  of  the  plasma  in 

order  to  correctly  simulate  MHD.  The  scalar  conductivity,  er,  is  defined  as: 

nee2 

a  =  - 

mevm 

where  ne  is  the  electron  number  density,  me  is  the  electron  mass,  and  e  is  the  Coulom- 
bic  charge  where  the  electron  collisional  frequency,  um  is  defined  as: 

/OO 

de  Qes(e)  v  •  / (e)  (2.28) 

^ - v- - - 

^en 
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where  Qei  and  uei  is  the  electron-ion  elastic  collisional  cross-section  and  collisional 
frequency,  Qes  and  ves  is  the  electron-neutral  elastic  collisional  cross-section  and 
collisional  frequency  of  the  sth  species,  ve  is  the  electron  thermal  velocity,  and  / 
is  the  electron  distribution  function.  Electron- neutral  cross-section  data  are  found 
in  [41]  while  the  electron-ion  cross-section,  commonly  referred  to  as  the  Coulombic 
cross-section  is  given  by  [42]: 


^  2.87  xlO"14,  . 

Qei  tTlA. 

I[eV] 


(2.29) 


where  the  Coulombic  logarithm  is  InA  =  13.57  +  1.5  log  XLyi  —  0.5  log  ne.  In  the 
regime  where  the  fluid  is  strongly  ionized,  a  >  10~4,  the  electron-ion  collision  term 
of  Equation  2.28  dominates,  and  the  conductivity  scales  as  follows: 


a  ~  Te5/2 


(2.30) 


2.3.3  Two-Temperature(2T)  Relaxation 

When  a  plasma  is  rapidly  heated  by  a  shock,  radiation,  or  other  process,  the  heavy 
particle  and  the  electron  temperatures  can  be  altered  from  equilibrium  and  must 
undergo  a  series  of  elastic  collisions  to  return  them  to  equilibrium,  or  thermalize.  If 
the  collisional  time  scale  is  significantly  faster  than  the  fluid  time  scale,  the  heavy 
particles  and  electrons  can  be  assumed  to  be  in  equilibrium.  But,  when  the  collisional 
time  scales  are  much  slower  than  that  of  the  fluid,  the  temperatures  can  be  treated 
completely  separately.  In  the  event  the  time  scale  associated  with  thermalization  is 
on  the  order  of  the  fluid  time  scale,  the  finite  rate  of  elastic  energy  relaxation  must 
be  taken  into  account.  The  relaxation  of  the  electron  energy  is  as  follows  [42]: 

d~w  =  Vm  (i t)  \nM  (Th  ~ T,)  (2-31) 
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where  Th  is  the  translational  temperature  of  the  heavy  particles.  Equation  2.31  can 
be  recast  in  terms  of  electron  entropy  as  follows: 


Se  = 


dSe 

dt 


=  Z4, 


(jr)  ln‘kb  (Th  ~  Te)  (7e  “ 1}  pl_7‘  (2-32) 


2.4  MHD  Transport 


In  order  to  properly  characterize  a  non-ideal  system  in  the  present  of  an  imposed 
magnetic  held,  one  must  account  for  both  the  convective  and  diffusive  transport  of 
the  magnetic  held.  Let  us  re-examine  the  formulation  of  the  temporal  evolution  of 
the  magnetic  held  expressed  in  Equation  2.10: 


dB 

~dt 


/i0(T 


-VxVxB  —  VxuxB 


- v - 

diffusive 


convective 


Starting  from  2.1  and  the  description  of  the  source  term  from  Eqn.  2.20,  one  can 
define  a  system  which  describes  the  evolution  of  the  magnetic  held  and  magnetic 
pressure  due  to  diffusive  transport  in  the  divergence  form  as  follows: 


9B 

dt 

dPB 

dt 


=  V- 
=  V- 


— VB 

MO  O' 

—  VT 

fi0a 


(2.33) 


where  Pb  is  the  magnetic  pressure  which  was  previously  defined  as  Pb  =  yyp  Now 
let  us  rewrite  Eqn.  2.33  in  a  hux  formulation  of  the  form, 


Qt  =  V  ■  F" 


(2.34) 
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where  F"  represents  the  magnetic  diffusive  flux, 


(  Bx  \ 

(  VHX  \ 

By 

F"  -  1 

VHy 

Bz 

fi0a 

VBZ 

\  Pb  ) 

<1 

(2.35) 


This  description  of  MHD  diffusive  transport  is  particularly  useful  when  the  con¬ 
ductivity  of  a  fluid  is  finite,  as  in  the  problems  represented  in  this  thesis.  The  limits 
of  ideal  MHD  and  perfectly  resistive  MHD,  where  cr  — »  oo  and  cr  «  0,  respectively 
will  also  be  investigated. 
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Species 

^  [feg-if] 

Cv  t  kg—K  ] 

e°  [fc5] 

h2-o2 

692.8 

2.425  x  103 

0 

h2o 

461.9 

2.079  x  103 

1.344  x  107 

Table  2.1:  List  of  thermodynamic  properties  of  a  stoichiometric  mixture  of  H2  —  02l 
reactant,  and  H20 ,  product,  for  a  simple,  single-step  reaction. 


Species 

Bi,  [K\ 

Cs 

45,141 

K 

50,364 

Na 

59,647 

Li 

62,548 

o2 

139,834 

h2o 

146,217 

0 

157,937 

co2 

167,105 

h2 

181,030 

n2 

181,030 

Ar 

182,887 

He 

285,239 

Table  2.2:  Listing  of  characteristic  temperatures  of  first  ionization  of  selected 
species  [43]. 
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1  /T initial  [\/K] 


Figure  2.1:  Computed  induction  delay  as  a  function  of  initial  temperature  at  various 
pressures  P  for  the  stoichiometric  H2-a.i1  reaction  with  complex  kinetics.  Curves  can 
be  fitted  as  t  ~  a(P)e^r1 . 
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CHAPTER  3 


Numerical  Methods 

The  following  chapter  will  list  and  describe  the  numerical  methods  utilized  to  solve 
the  governing  equations  described  in  the  previous  section. 

3.1  Operator  Splitting 

As  shown  in  Equation  2.4,  it  is  often  necessary  to  solve  the  LHS  with  one  solver  while 
solving  the  various  RHS  source  terms  with  other  solvers.  For  example,  one  might 
wish  to  solve  a  problem  which  involves  convection,  LHS,  as  well  as  chemical  kinetics, 
RHS,  and  utilize  different  solvers,  or  operators ,  Cconv  and  JCchemi  respectively.  By 
utilizing  these  operators,  the  solution  of  Q  at  t  =  tn+1,  where  tn+1  —  t0  +  At,  can  be 
determined  as  follows: 

AQ chem  =  CchemQn  -  Q" 

AQconv  =  £convQn  -  Qn  (3.1) 

Qn+1  =  Q”  +  AQchem  +  AQconv 

Figure  3.1  illustrates  the  contribution  of  the  two  operators  in  attaining  the  solution 
at  t  —  tn+1. 

When  stronger  coupling  is  needed  amongst  multiple  operators,  either  for  stability 

or  accuracy,  Strang  splitting [44]  can  be  employed.  Strang  splitting  the  convective 
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and  chemical  kinetic  operator  of  the  previous  example  would  go  as  follows: 


Q 

Qn+1 


=  Cchem  Q 
=  Cconv  Q 


n 


which  can  be  simply  expressed  in  the  compact  notation: 


Qn+1  _  £conv £chemQn 


(3.2) 


Figure  3.2(a)  illustrates  how  the  chemical  kinetic  solver  operators  on  the  solution  at 
t  =  tn,  and  the  convective  solver  operates  on  the  updated  solution  in  order  to  obtain 
the  solution  at  t  =  tn+1.  Since  Strang  splitting  is  not  a  commutative  operation,  the 
ordering  of  the  operators  can  indeed  yield  unique  solutions, 


£chem £C.onv qti  -J-  £C.onv £C/iemQ 


conv  nchem / 


(3.3) 


Figures  3.2(a)  and  3.2(b)  illustrate  how  commuting  the  operators  could  yield  a  dif¬ 
ferent  solution. 


3.2  Time  Step  Restrictions 

3.2.1  Convection 

For  linear  or  linearized  hyperbolic  equations,  the  Courant-Friedrich-Levy(CFL)  con¬ 
dition  ensures  that  information  does  not  propagate  further  than  one  grid  cell  in  a 
given  time  step: 

7/ A  r 

A tconv  <  (3.4) 

where  A tconv  represents  the  convection-limited  time  step,  v  is  the  Courant  number, 

for  linear  stability  0  <  v  <  1,  Ax  is  the  grid  cell  size,  A  is  the  wave  speed  of  the 
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convected  information.  For  a  purely  convective  simulation,  where  there  is  no  kinetic 
or  diffusive  processes,  the  CFL  would  be  sufficient,  for  stability  of  the  time  integration. 

3.2.2  Kinetics  Resolution 

Solving  kinetics  presents  a  different  challenge  in  that  the  reaction  time  scale  must  be 
taken  into  consideration  in  order  to  determine  the  proper  time  step  restriction.  For 
the  simple  case  of  an  explicit  single  step  reaction  with  the  form,  6jp  =  pfKe~9^T ,  the 
time  step  restriction  would  be  determined  as, 

A tchem  <  \K~1e9/T\  (3.5) 

In  the  case  of  an  explicit  complex  reaction  of  the  form  of  Equation  3.41, 

A  tchem  <  min  f  (3.6) 

\UJk  oje  J 

Since  the  kinetic  source  term  of  Equation  2.3  is  often  extremely  stiff  and  A tchem  ~  0 
when  pk  ~  0,  an  implicit  approach  will  be  taken.  When  a  strictly  implicit  numerical 
scheme  is  employed,  the  time  step  restrictions  for  stability  of  Equation  3.6  are  no 
longer  required.  But  even  though  restrictions  for  stability  are  no  longer  required, 
there  still  exist  accuracy  considerations.  The  following  accuracy  time  step  restriction 
was  introduced  in  the  present  study: 

A  tchem  =  min  ^ecj£^,  6t  |^j)  (3-7) 

Where  p*k  is  the  greater  value  of  pk  and  a  floor  density,  ec  is  the  maximum  species 
creation/destruction  fraction,  and  is  the  maximum  temperature  cooling/heating 
fraction. 
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3.2.3  Diffusion 


When  solving  diffusive  transport  problems,  which  are  expressed  as  parabolic  partial 
differential  equations  of  the  form,  ut  =  ~^uxx,  the  explicit  time  step  restriction  is 
quite  similar  to  that  of  the  convective  transport.  The  explicit  time  step  restriction 
is  given  by, 


0.5Aa;2 


(3.8) 


where  /i  represents  the  “kinematic  viscosity”.  This  time  step  can  become  quite 
constraining  with  the  increase  in  the  diffusivity,  /x,  or  with  moderate  grid  spacing 
reduction  (At  ~  Ax2).  For  these  reasons,  spatially  implicit  numerical  schemes  will 
be  utilized  to  solve  problems  of  this  type. 

3.3  Explicit  Runge  Kutta  Scheme 

In  order  to  obtain  stable  high  order  convergent  solutions  in  time,  a  TVD  Runge-Kutta 
(RK)  time  integrator  is  used.  This  particular  version  of  the  Runge  Kutta  family  was 
implemented  by  Shu  and  Osher  in  [45].  For  3rd  order  Runge-Kutta  (RK3),  V  j 


(3.9) 


where  Lj  represents  the  flux  into  the  jth  grid  cell  and  represents  the  conserved 
variables  of  the  kth  step  of  the  Runge-Kutta  integration.  Equation  3.9  will  serve  as 


the  fluid  convection  operator,  Cconv ,  for  the  remainder  of  this  study,  unless  otherwise 
stated,  while  the  temporally  and  spatially  integrated  flux,  Lj ,  will  be  determined  by 
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Riemann  solvers  which  will  be  discussed  in  the  next  sections. 


3.4  Hyperbolic  Solvers  —  Approximate  Riemann  Solvers 


In  the  present  study,  it  is  imperative  that  the  fluid  convection  solver  be  able  to 
capture  shocks  without  introducing  dispersion  or  excess  dissipation  into  the  solution. 
Let  us  now  recast  the  LHS  of  our  governing  equation,  Equation  2.1,  into  its  integral 
or  conservative  form: 


d_ 

di 


Q  dV  + 


F  dS  =  0 


(3.10) 


We  will  next  assume  a  given  grid  cell  is  uniform  in  its  properties,  and  can  be  ap¬ 
proximated  with  its  cell-averaged  value.  Then  Equation  3.10  can  take  the  following 
form: 

f  4LFA=0  t3-11) 

s 

where  Fs  and  dAs  are  the  areas  and  normal  fluxes  of  the  sth  face  of  a  given  grid  cell. 
As  a  building  block,  we  will  utilize  Roe’s  scheme [46]  to  solve  for  the  flux  at  the  face, 
Fj_|_i/2  as  follows: 


Fh-1/2 


-CFr 
2  v 


-RAL  (Q# 


Q  i 


(3.12) 


where  F  and  Q rl  represent  the  normal  fluxes  and  conserved  variables  of  the  right 
and  left  side  of  the  face,  respectively,  L  is  the  matrix  of  eigenvector,  R  is  the  inverse 
of  the  matrix  of  eigenvector,  and  A  are  the  HLLE[47]  conditioned  eigenvalues.  The 
eigen-system  is  discussed  in  more  detail  in  Appendix  B.  For  the  1st  order  spatial 
accuracy  convergent  Roe  scheme,  the  right  fluxes  and  conserved  variables  are,  F;+1 
and  Qi+i,  respectively  and  the  left  fluxes  and  conserved  variables  are,  F,  and  Q,, 
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respectively.  In  the  proceeding  sections,  the  numerical  schemes  used  to  calculate 
high  order  spatially  convergent  conserved  variable  solutions  will  be  discussed. 


In  order  to  achieve  a  spatially  high  order  convergent  interpolation  of  the  conserved 
variables  at  the  cell  interface,  Qr,l,  the  governing  system  of  equations,  Equation  2.1, 
must  be  linearized  in  the  following  manner: 


Qi  +  Fn  —  0 

Q  t  +  AQn  =  0 

LQt  +  LRALQn  =  0 

Wt  +  AW„  =  0 


(3.13) 


where  ()n  represents  the  spatial  derivative  in  an  arbitrary  direction,  the  convective 
flux  jacobian  is  A  =  and  the  characteristic  variable  array,  W  =  ( W\,W2 ,  ■■■)T  is 
defined  as  the  projection  of  the  conserved  variables,  W  =  LQ,  and  by  definition 
LR  =  RL  =  I.  Now  that  the  governing  equations  have  been  linearized  with  the  pro¬ 
cess  shown  in  Equation  3.13,  it  can  now  be  expressed  as  a  system  of  scalar  hyperbolic 
differential  equations: 

wt  +  A  wn  =  0  (3-14) 


where  the  eigenvalues  are  the  diagonal  components  of  the  matrix  of  eigenvalues, 
A i  =  A jtj.  After  using  one  of  the  high  order  spatially  convergent  methods,  which 
will  be  discussed  in  greater  detail  in  the  preceding  sections,  to  approximate  the 
characteristic  variable  solution  at  the  cell  interface,  W the  characteristic  variables 
can  be  projected  back  to  its  component  form  using  the  following  operation: 


Q  r,l  —  RW  r}l  (3.15) 

The  updated  conserved  variables  determined  from  Equation  3.15  are  then  used  to 

calculate  the  interface  flux,  FRL,  as  well  as  construct  a  new  eigensystem. 
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3.4.1  Weighted  Essentially  Non-Oscillatory  (WENO)  Schemes 

In  Weighted  Essentially  Non-Oscillatory  (WENO)  schemes,  first  introduced  by  Liu, 
Osher,  &  Chan  [48] ,  high  spatial  order  of  convergence  is  achieved  where  the  solution 
is  smooth  and  a  spatial  convergence  of  no  greater  than  0(1)  near  a  discontinuity 
in  the  solution.  In  the  present  study,  we  utilized  a  5th  order  spatially  convergence 
variant  of  WENO  which  weights  the  contribution  of  three  stencils,  illustrated  in 
Figure  3.3.  Here  we  will  describe  how  u>l, i+1/2  is  computed  in  [48]  on  the  basis  of  the 
ENO  stencil  [49].  For  simplicity,  the  “L”  subscript  will  be  dropped.  The  formula  for 
the  right  characteristics  are  symmetric  and  will  only  be  shown  when  they  vary  from 
the  left  characteristics. 

The  rth  order  ENO  scheme  chooses  the  “smoothest”  stencil  from  r  candidate 
stencils  to  approximate  Wi+ 1/2.  In  the  case  of  r  =  3,  the  stencil  Sk,  where  k  E  [0,  2], 
shown  in  Figure  3.3,  happens  to  be  chosen  as  the  ENO  interpolation  stencil,  the 
W'-order  ENO  approximation  of  wy+1/2  to  produce 

Wk  Q 'k  i^i+k— r+1  j  •••;  ^ i+k )  (3.16) 

where 

r—  1 

qrk(9o,--,gr-i)  =  ^ark>lgi  (3.17) 

1=0 

Here  ark  z,  0  <  k,  l  <  r  —  1,  are  constant  coefficients,  which  are  provided  in  Table  3.1. 
Using  the  smoothest  of  the  rth  stencil  would  be  desirable  near  a  discontinuity,  but  in 
smooth  regions,  information  from  all  stencils  can  be  used  in  the  final  solution.  Thus, 
in  smooth  regions,  it  would  be  desirable  to  combine  the  ENO  stencils  in  a  manner 
that  will  generate  a  higher  than  r  order  solution.  As  shown  in  [48],  one  can  use  all 
of  the  r  candidate  stencils,  which  all  together  would  contain  (2r  —  1)  grid  values  of 
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wi+ 1/2  to  give  a  (2 r  —  l)th-order  approximation  of  w 

Wk  =  q^Si1  (Wi-r+l,  Wi+r-l)  (3.18) 

which  is  the  solution  of  a  (2r  —  l)th-order  upstream  central  scheme.  Since  high 
order  upstream  central  schemes  (in  space),  combined  with  high  order  Runge-Kutta 
methods  (in  time)  are  stable  and  dissipative  under  appropriate  CFL  numbers,  they 
are  convergent.  Using  this  fact,  one  can  now  use  the  (2r  —  l)th-order  upstream 
central  scheme  in  smooth  regions  and  the  rth  order  ENO  scheme  near  discontinuities. 
As  shown  in  Equation  3.16,  each  of  the  stencils  can  approximate  wi+i/2.  If  the  stencil 
is  smooth,  an  rth  order  approximation  of  the  stencil  can  be  recovered,  but  if  the  stencil 
is  discontinuous,  a  less  accurate  or  inaccurate  approximation  would  be  recovered.  So 
WENO  assigns  a  weight,  u A,  to  each  of  the  candidate  stencil  Sk,  where  k  E  [0,  r  —  1], 
and  uses  these  weights  to  combine  the  r  different  approximations  to  obtain  the  final 
approximation  of  the  solution  as: 

r— 1 

^0+1/2  ^  (  0 kQk  i^i+k— r+li  (3.19) 

k= 0 

where  qrk  is  defined  in  Equation  3.17.  To  achieve  essentially  non-oscillatory  properties, 
WENO  requires  that  the  weights  adapt  to  the  relative  smoothness  of  w.  Discontinu¬ 
ous  stencils  contributions  should  be  assigned  weights  of  zero.  In  the  smooth  regions 
the  weights  should  be  adjusted  so  the  upstream  central  scheme,  Equation  3.18,  is 
recovered. 

In  the  present  study,  the  5th  order  WEN0(WEN05)  scheme  (r  =  3)  is  utilized 
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as  prescribed  by  Jiang  and  Shu  in  [45].  The  stencils,  Sk,  are  calculated  as: 

Soiwi-2,  Wi^Wi)  =  |  (2wi_2  -  7wi-i  +  11  Wi) 

Siiwi-^w^Wi+i)  =  1  +  5wi  +  2wi+1)  (3.20) 

S2(wi,  wi+uwi+2)  =  |  (2 Wi  +  5wi+1  -  lwi+2) 

And  the  right  solution  is  stencils,  S*  are  calculated  as: 

So(wi+ 3,wi+2,wi+1)  =  1  ( 2wi+3  -  7wi+2  +  llwi+1) 
SP(wi+2,Wi+1,Wi)  =  1  {~lwi+2  +  5wi+1  +  2wi)  (3.21) 

S?(wi+i,Wi,  Wi_ i)  =  1  (2wi+i  +  5 Wi  -  ltWi-i) 

The  smoothness  of  each  stencil  is  then  calculated  as: 

I S0  =  (wi-2  ~  2w*_i  +  Wi  f  +  1  (wi_ 2  ~  4wj_i  +  3 w^2 

I  Si  =  §  {wi- 1  -  2  Wi  +  wi+1  f  +  1  («;*_!  -  wi+i)2  (3.22) 

(Wi  -  2wi+i  +  tty+2)2  +  \  (3 Wi  -  Awi+i  +  wi+2)2 
The  smoothness  of  each  stencil  for  the  right  solution  is  then  calculated  as: 

IS*  =  If  (wj+3  -  2wi+2  +  wi+i )2  +  1  (uy+3  -  4wi+2  +  3wi+i)2 

I  Si  =  If  (wy+2  -  2rci+i  +  wy)2  +  \  ( wi+2  -  in*)2  (3.23) 

IS*  =  If  (wi+i  -  2 Wi  +  Wi-if  +  1  (3iyi+i  -  Awi  +  wy-i)2 

We  can  now  calculate  the  new  weights,  £>[,  based  on  the  smoothness,  ISk,  and 
the  optimal  weights,  ,  which  are  defined  as  uopt  =  (^,  A)  for  approximation 
of  w  by  using  the  procedure  from  [45] , 


u 


opt 


OJu  = 


(3.24) 


k  (e  +  ISk )2 

where  e  =  ICE6  is  there  to  guarantee  non-singular  behavior.  These  new  weights,  w'k, 

are  then  normalized  to  become  the  final  WENO  weights: 
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(3.25) 


where  ^2ku)k  =  1-  The  weighted  solutions  from  each  stencil  are  then  summed,  in  the 
same  manner  as  Equation  3.18,  to  form  our  WEN05  approximated  solution: 


2 


Wi+i/2  =  “kvl 


(3.26) 


k= o 


After  at  characteristics  at  all  of  the  cell  interfaces  are  calculated,  the  charac¬ 
teristics  are  then  projected  back  to  real  space  using  Equation  3.15  then  the  final 
flux  at  the  cell  interface,  F i+1/2,  is  determined  by  Equation  3.12  using  the  updated 
approximations  to  the  conserved  variables. 

3.4.2  Monotonicity  Preserving  (MP)  Schemes 

The  Monotonicity  Preserving  (MP)  scheme,  first  introduced  by  Suresh  &  Huynh[50], 
uses  a  high  order  spatially  convergent  reconstruction  of  the  interface,  the  original 
value,  then  limits  this  solution  in  order  to  obtain  the  final  interface  value.  We  will 
adopt  the  notation  of  the  previous  section  and  drop  the  “L”  subscript  for  solutions 
at  the  face.  For  the  5th  order  MP  scheme  (MP5),  which  will  be  used  throughout  the 
present  study,  the  original  value,  w^/2,  is  given  as: 


w?+i/2  —  (2wi_2  -  13uy_i  +  47 Wj  +  27-uy+i  -  3wi+2)  /60  (3.27) 


To  find  the  MP5  solution,  a  few  constraints  must  be  satisfied.  The  first  constraint 
is  monotonicity- preservation,  which  Suresh  &  Huynh  [50]  define  as  the  upper  limit, 


wUL  =  Wi  +  a(wi  -  Wi- 1) 


(3.28) 
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where  a  =  2,  represents  a  parabolic  reconstruction.  The  second  constraint,  accuracy- 
preserving,  is  accomplished  by  bounding  the  solution  with  the  median  and  large  cur¬ 
vature  solutions.  The  median  (denoted  by  superscript  “MD”)  value  at  the  interface, 


Xi+1/2,  is  given  by 

mMD  =  WAV  _  1  ^  (3.29) 

While  the  large  curvature  (denoted  by  superscript  “LC”)  value  at  the  interface  is 
given  by 

wLC  =Wi  +  ^(wi-  Wi- 1)  +  i/2  (3.30) 

where  (3  —  4,  wAl  is  the  average  solution,  di  is  the  curvature,  and  df^1j2  is  the 
minmod  approximation  of  the  curvature  at  the  zone  boundary,  all  of  which  are 
defined  as  follows: 

WAV  =\(Wi+Wi+i) 

di  =  i  +  wl+x  -  2 Wi  (3.31) 

diif/2  =  minmod  (di,  di+1) 

The  superscript  “MM”  indicates  the  use  of  a  minmod  function.  Suresh  &  Huynh [50] 
recommended  the  use  of  a  slightly  more  restrictive  curvature  measure  than  df^1 
which  is  given  by: 


dfii/2  =  minmod  (Adi  ~  di+1,  Adi+l  -  di,  dh  di+1)  (3.32) 

For  the  MP5  scheme  d^/2  =  di+i/2  =  dtl\/2-  Now  that  the  mechanisms  for  the  two 
constraints,  monotonicity- preserving  and  accuracy- preserving,  have  been  stated,  the 
minimum  and  maximum  value  of  the  solution,  wMIN  and  wMAX ,  respectively,  are 
given  by: 

wmin  _  max  ^min(ryj,  min(iUj+1,  wMD)),  min(ryj,  m\n(wUL ,  wLC))^ 

(3.33) 

wMAX  =  min  [max(tCj,  max(tCj+1,  wMD)),  max('u,'„  ma x(wUL,wLC))] 
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The  solution  at  the  face,  wi+i/2,  can  now  simply  be  expressed  as: 

wi+ 1/2  =  median  wMIN ,  wMAX)  (3.34) 

The  5th  order  convergent  Monotonicity  Preserving  WENO  (MW5)  scheme  by 
Balsara  &  Shu  [51],  which  is  based  on  the  MP5  scheme,  has  the  smooth  solution 
of  WEN05,  but  is  strictly  monotonicity-preserving  near  discontinuities.  Balsara  & 
Shu  were  able  to  demonstrate  that  the  solutions  to  ideal  MHD  simulations  were  far 
superior  to  that  of  MP5  and  WEN05.  The  form  of  the  MW5  solution  is  the  same 
MP5,  Equation  3.34,  but  the  key  difference  is  in  how  MW5  is  used  to  determine  the 
original  value,  wOR,  and  the  curvature  at  the  median  and  large  curvature  values, 
dAID  and  dLC ,  respectively.  By  using  the  WENO  stencils  for  r  =  3,  Equation  3.20 
with  WENO  coefficients  from  Table  3.1,  and  the  optimal  weights,  the  MW5  original 
value,  wOR,  is  recovered: 

w°h  =  (6wy_2  —  27tCj_i  +  65 wy  +  YJwi+i  —  wi+ 2)  / 60  (3.35) 

In  an  attempt  to  filter  out  extremal  features  with  small  domains  of  support,  while 
keeping  extremal  features  with  large  domains  of  support  intact,  =  dR^,2  = 

^i+i/2>  where 

d^X/2  =  minmod  (Adt  -  di+1,  Adi+1  -  dt,  di,  di+1,  dj_i,  di+2)  (3.36) 

3.4.3  Advection-DifFusion-Reaction  (ADER)  Schemes 

The  Advection-Diffusion-Reaction  (ADER)  schemes  of  Titarev  &  Toro  [52]  utilize 
high  order  spatial  derivatives  calculated  by  the  underlying  scheme  to  generate  the 
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temporal  derivatives  using  the  Cauchy-Kovalevskaya  procedure: 


dtw  =  —A  dxw 
dttw  =  A  2dxxw 

dtttw  =  -X3dxxxw  (3.37) 


d^w  =  (— A)fc<9 


Where  A  is  the  characteristic  wave  speed.  With  the  high  order  temporal  derivatives 
generated  from  Equation  3.37,  a  simple  Taylor  series  expansion  is  performed  to  ac¬ 
quire  a  higher  order  temporally  and  spatially  convergent  scheme.  But  first,  we  must 
take  the  temporal  series  expansion  at  the  interface, 

m—  1  fc 

Wi+l/2  =  W  +  ^  (3-38) 


Where  r  is  the  time  step  size  and  w  =  u>(oy+i/2,  0+)  is  the  approximation  of  the 
solution  at  the  interface  from  the  underlying  scheme.  In  the  present  study,  we  wish 
to  achieve  3rd  order  convergence  in  time,  so  we  will  start  from  Equation  3.38  with 
m  =  3  to  give, 


Wi+l/2 


w  +  dtWT  +  duw 


2 


(3.39) 


By  performing  the  Cauchy-Kovalevskaya  procedure  on  Equation  3.39,  the  solution, 
Wi+ 1/2,  will  become  a  function  of  the  time  step  size  and  spatial  derivatives  of  the 
approximate  solution  as  follows: 


U>i+ 1/2 


W 


XdxuiT  +  X2dxxw— 


(3.40) 


Equation  3.40  will  serve  as  the  general  form  of  the  temporally  3rd  order  convergent 
ADER  (ADER3)  scheme.  Since  the  expensive  Runge-Kutta  time  integration  steps 

56 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


are  no  longer  required,  ADER  schemes  are  extremely  efficient  and  well  suited  for 
parallel  computation. 

By  combining  the  5th  order  spatially  convergent  WEN05’s  reconstruction  to  cal¬ 
culate  the  approximate  solution  and  its  spatial  derivatives,  w  and  dxiv,  respectively, 
with  the  3rd  order  temporally  convergent  ADER3,  we  form  the  3rd  order  temporal 
and  5th  order  spatial  convergent  ADER-WENO  (AW5)  scheme  [53].  At  discontinu¬ 
ities,  di^w  =  0  V  k  >  1  is  satisfied  in  order  to  prevent  spurious  oscillations. 


3.5  Point-Implicit  Euler 


The  finite  rate  kinetic  systems  in  the  present  study  are  extremely  stiff.  When  solving 
a  stiff  ODE,  it  is  often  beneficial  to  solve  the  problem  implicitly.  The  stiff  chemical 
kinetics  ODEs  will  be  expressed  in  the  form  of  Equation  2.3: 


dQ 

dt 


kinetic 


(3.41) 


where  Q  and  S2 kinetic  are  defined  by  Equation  2.16.  The  1st  order  point-implicit  Euler 
will  be  utilized  to  solve  problems  of  this  particular  form.  First,  let  us  discretize  the 
time  into  uniform  intervals  of  size  At  and  denote  tn  =  t0  and  tn+1  =  t0  +  At.  Upon 
discretization,  a  Taylor  series  expansion  can  be  performed  on  Equation  3.41  and  is 
expressed  as  follows: 


A  Q 
At 

A  Q 
At 

A  Q 
At 

A  Q 
At 


=  nn  +  f  rf  a  t 

=  Vn  +  ijViQAt 

=  +  %  AQ 
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Then  we  shall  proceed  and  solve  for  the  change  in  the  conserved  variables,  A Q,  and 
arrive  with  the  following  form: 


AQ  =  I 


on 

~dQ 


-i 


■At  ft  At 


(3.43) 


From  here  we  arrive  at  the  solution  to  our  implicit  formulation  as: 


Qn+1  =  Qn 


T  <90  A  ,  .  n  A 

a  At 


(3.44) 


This  can  also  be  expressed  in  the  operator  form:  Qn+1  =  CchemQn,  where  the  oper¬ 
ator  is  defined  as: 

I  an"  \ _1  ft  n” 

(3.45) 


,  I  dn  \  dQ, 

Cchem  =  1  +  j  _  __At  }  ——At 


dQ 


dQ 


3.6  Spatial-Implicit  Euler 


Diffusive  MHD  transport  will  often  have  much  stricter  explicit  time  step  restrictions, 
Equation  3.8,  than  that  of  convective  transport,  Equation  3.4.  Since  the  maximum 
allowable  explicit  diffusive  time  step  is  determined  by  At  <  0.5Ax2aminfio,  it  becomes 
apparent  that  as  a  — »  0  then  At  — >  0.  This  can  become  quite  cost  prohibitive,  so  in 
order  to  ensure  stability  with  a  non  prohibitive  time  step,  an  implicit  time  marching 
scheme  is  utilized. 

Before  we  can  cast  the  diffusive  MHD  transport  into  an  implicit  formulation, 
we  must  start  from  the  magnetic  diffusion  flux  formulation  of  the  RHS,  Equation 
2.34,  and  define  the  magnetic  field  diffusion  Jacobian,  Av ,  for  the  dimensionally  split 
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magnetic  diffusion  as 

9Q  =  ^/  3Q\ 

dt  dx\  x  dx  ) 

dFu 

where  the  magnetic  field  diffusion  Jacobian,  Ax  = 


A 


V 

X 


( 


1 

/i0cr 


1 

0 

0 


0  0  (0 
1  0  0 
0  1  0 
0  0  1  J 


(3.46) 


Assuming  a  one-dimensional  discretization  on  a  uniformly-spaced  grid,  the  spatial 
derivatives  can  be  approximated  by  finite-differences  and  the  subscript  in  Ax  will  be 
ignored,  thereby  reducing  the  PDE  to  a  system  of  ODE’s, 


dQj 

dt 


Ax2 


(3.47) 


the  system  of  equations  can  be  written  in  matrix  notation  as 


dQ 

dt 


=  <lQ 


(3.48) 


where  $  is  a  tridiagonal  matrix  and  Q  is  the  spatial  vector  of  the  conserved  element 
array, 


/ 


$ 


1 

Ax2 


V 


\ 


/  :  \ 


Qi— 1 

Q, 

Qt+i 


/  \  ]  ) 

(3.49) 
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Now,  we  can  determine  the  implicit  formulation  by  evaluating  the  RHS  of  Eqn. 
3.48  at  time  t  +  At, 

Q"+1  =  Qn  +  At  $  Q"+1  (3.50) 

By  applying  backwards  Euler  to  Equation  3.50,  one  can  recover  the  operator  form, 
Qn+1  =  CdJff  Qn,  where  the  MHD  diffusion  matrix  operator  for  the  Q-direction  is 

Cfff  =  (l-  (3.51) 

In  order  to  apply  the  diffusive  MHD  transport  in  all  directions,  one  might  perform 
the  following  operation: 

Qn+1  =  CdiffCdiffCfff  Qn  (3.52) 

Upon  inspection  of  Equation  3.49,  £*1/  js  merely  the  inverse  of  a  tridiagonal 
system  of  equations.  Rather  than  applying  a  scheme  with  a  relatively  high  computa¬ 
tion  cost,  i.e.  GMRes  or  Gaussian  Elimination,  one  can  exploit  the  fact  the  system 
is  tridiagonal  and  implement  Thomas’  Algorithm,  which  is  discussed  in  more  detail 
in  Appendix  F.l. 

The  diffusive  MHD  operators,  Cfx^ ,  use  line  relaxation  to  proceed  in  time.  If 
the  grid  and  fluid  properties,  i.e.  conductivity,  grid  resolution,  etc.,  are  not  spa¬ 
tially  uniform  the  diffusive  MHD  transport  will  have  different  time  scales  for  each 
directional  sweep.  When  this  occurs,  it  becomes  necessary  to  split  a  given  operator, 
jCx,  in  time  using  the  Strang  operator  splitting  technique  demonstrated  in  Section 
3.1.  The  superscript  ‘diff’  has  been  discarded  for  the  remainder  of  this  section.  One 
permutation  of  applying  multiple  spatial  operators  to  a  generic  2D  diffusive  MHD 
problem  would  go  as  follows: 

Qn+1  =  £f /2  £f  £f /2Qn  (3.53) 
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which  is  illustrated  in  Figure  3.4(a).  In  order  to  avoid  developing  a  bias  toward  a 
particular  direction,  after  applying  Equation  3.53  one  should  use  the  following: 

gn+2  =  £At/2  £A t  £Ai/2gn+l  (3.54) 

which  is  illustrated  in  Figure  3.4(b).  A  similar  permutation  can  be  performed  when 
solving  a  3-D  diffusive  MHD  problem. 
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/  =  0 

l  =  1 

/  =  2 
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3/2 
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1/3 

-7/6 

11/6 

1 

-1/6 

5/6 

1/3 

2 

1/3 

5/6 

-1/3 

Table  3.1:  WENO  Coefficients,  ark  l 


Q 


Figure  3.1:  Operator  splitting  of  two  generic  operators,  Cchem  &  £com'. 
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Figure  3.2:  Demonstration  that  with  Strang  splitting  coupled  operators,  the  ordering 
of  the  operators  can  significantly  change  the  solution. 
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Figure  3.3:  WENO  stencil  for  r  =  3. 
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Figure  3.4:  Strang  operator  splitting  in  time  for  three  operations  illustrated  with 
different  ordering  of  the  operators,  Cx  and  Cy. 
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CHAPTER  4 


Verification  of  Numerical  Schemes 


The  numerical  schemes  used  in  the  present  study  must  be  validated  and  verified  to 
ensure  there  accuracy  and  proper  usage.  The  1-dimensional,  2-dimensional,  and  3- 
dimensional  forms  of  the  hydrodynamic  and  magnetohydrodynamic  governing  equa¬ 
tions  assume  the  form  of  Equation  2.1: 


DQ 

dt 


+  Vn  •  F 


n 


For  the  1-D  subset  of  the  governing  equations,  the  variation  in  the  y-  and  z- 
dimensions  are  set  to  zero,  ^  =  0  and  =  0,  respectively.  While  in  the  2-D  subset 
of  the  governing  equations,  there  is  no  variation  in  the  z-direction,  =  0. 


4.1  Inviscid  Hydrodynamics 

The  hydrodynamic  test  cases  will  show  that  the  physical  wave  (acoustic  and  entropy) 
speeds  present  in  the  simulation  are  properly  captures  and  that  the  approximate 
Riemann  solution  matches  the  proper  jump  conditions  in  the  event  of  a  shock.  In 
the  hydrodynamic  limit  the  conserved  variables,  Q,  and  flux,  F,  are  expressed  as 


65 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


Equation  2.5: 


Q  = 


p  N 

pu 

\  E  ) 


F  = 


\ 


pUn 

pu  ■  un  +  Pn 
y  (E  +  P)  un  J 


4.1.1  Sod’s  Shock  Tube 


There  are  a  few  classic  test  cases  performed  utilizing  our  Riemann  solvers,  in  order 
to  test  their  ability  to  resolve  shocks,  contact  discontinuities,  and  rarefactions  as  well 
as  the  interactions  of  these  structures.  The  first  test  case  was  introduced  by  Sod  [54], 
known  presently  as  Sod’s  shock  tube.  We  will  adopt  the  same  initial  and  boundary 
conditions.  The  left  and  right  states  are 


(. Pl,uXjL,Pl )  =  (lkg/m3,0m/s,l  Pa)  x<0 

C Pr,ux,r,Pr )  =  (0.125  kg/m3,  0  m/s,  0.1  Pa)  x>0 


(4.1) 


where  the  adiabatic  index,  7  =  1.4,  unless  otherwise  stated,  and  over  a  domain 
x  G  (—5,5 ]mm.  Figure  4.1(a)  shows  the  solution  for  the  schemes  at  t  =  200  ps 
which  contains  a  single  shock,  contact  discontinuity  and  rarefaction.  The  figure 
clearly  shows  that  all  of  the  schemes  used  resolved  the  problem  reasonably  well 
compared  to  the  exact  solution.  Figure  4.1(b)  illustrates  MW5  solution  converging 
to  the  exact  solution  at  the  contact  discontinuity  as  Ax  is  decreased. 


4.1.2  Hydrodynamical  Interacting  Blast  Wave  Problem 

Blast  waves  are  generally  described  as  strong  and  rapid  release  of  energy  which  are 
often  characterized  by  regions  containing  drastic  temperatures  and  pressure  rises. 
For  the  second  test  of  our  ID  hydrodynamic  test  problems,  we  ran  the  interacting 
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blast  wave  problem  which  was  first  proposed  by  Woodward  and  Collcla  [55].  The 
problem  consist  of  the  left,  middle,  and  right  initial  states,  L,  M,  and  R  respectively, 
which  are  as  follows: 


(pL, 

^ x,l ?  Pl) 

=  (1,0, 103) 

x  <  0.1 

(, PM-,  Ux,Mi  Pm) 

o 

o 

1 

to 

0.1  <  x  <  0.8 

(. PRj 

ux,R,  Pr) 

=  (1,0, 102) 

x  >  0.8 

(4.2) 


where  x  G  (0,1].  Figure  4.2  shows  the  density  profiles  at  t  =  38  seconds  using 
the  MP5,  AW5,  and  MW5  schemes  with  greater  detail  in  Figure  4.3.  While  all 
of  the  schemes  resolve  the  shocks  at  x  ~  0.65  &  0.87  reasonably  well,  the  contact 
discontinuities  at  x  ~  0.75  &  0.8  are  slightly  better  resolved  by  MP5.  None  of  the 
schemes  used  artificial  compression  methods. 


4.1.3  Shock-Entropy  Wave  Interaction 

In  the  last  of  our  1-D  hydrodynamic  test  cases,  we  wish  to  test  our  numerical  schemes 
ability  to  resolve  smooth  flow  disturbances  which  is  of  particular  interest  because 
of  the  nature  of  the  problems  in  the  present  study.  The  Shu-Osher  problem  [56] 
has  been  extensively  used  to  simulate  a  Mach  3  shock  wave  interacting  with  an 
oscillatory  density  disturbance  which  generates  a  flow  field  with  a  combination  of 
smooth  structures  and  discontinuities.  The  initial  conditions  are  given  for  the  left 
and  right  state  as  follows: 

(, Pl,uX}L ,  Pl)  =  (3.857143,  2.629369, 10.33333)  x  <  0.8  ^ 

(. PRiUx,RiPr )  =  (1  +  0.2sm(57rx),0, 1)  x  >  0.8 

where  x  G  (—1, 1  ]m.  From  the  density  profile  in  Figure  4.4(a)  at  t  —  360 ms,  one 

can  see  that  three  schemes  resolved  the  entropy  disturbances  quite  well.  Upon  closer 
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inspection,  Figure  4.4(b)  clearly  shows  that  while  AW5  was  slightly  diffused  and 
MW5  slightly  amplified,  MP5  was  clearly  able  to  resolve  the  transient  entropy  waves 
the  best. 

4.1.4  Shock  Diffraction  Down  a  Backward  Facing  Step 

The  next  test  problem  describes  the  diffraction  of  a  Mach  2.4  shock  down  a  backward 
facing  step [57].  The  strong  rarefaction  generated  by  the  diffraction  at  the  90°  corner 
often  results  in  numerical  errors  described  by  over-expansion  and  negative  pressure 
for  many  Riemann  solvers  [58].  The  problem  is  simulated  using  a  resolution  of  Ax  = 
Ay  =  with  the  MW5  solver.  The  numerical  simulation  is  shown  side-by-side 
with  the  experimental  images  in  Figure  4.5.  The  numerical  solution  is  presented 
using  a  Schlieren-type  plot  as  prescribed  by  [58]  which  uses  density  gradients  in  an 
analogous  way  to  index  of  refraction  gradients,  which  makes  it  ideal  for  comparison 
with  experimental  images.  The  figure  shows  that  the  MW5  scheme  was  able  correctly 
reproduce  the  flow  features  in  the  region  of  the  rarefaction. 

4.1.5  Rayleigh- Taylor  Hydrodynamic  Instability 

In  the  next  test,  a  heavy  fluid  is  supported  by  a  lighter  fluid  in  a  gravitational  held, 
or  equivalent,  which  accelerates  the  heavier  fluid  into  a  lighter  fluid.  This  condition 
is  unstable  once  the  interface  between  the  two  fluids  is  perturbed.  The  instability  is 
known  as  the  Rayleigh- Taylor (RT)  instability.  Earlier  analytical  investigations  date 
back  to  the  detailed  analysis  given  by  Chandrasekhar [59]. 

In  the  initial  configuration,  two  fluids  with  a  prescribed  density  ratio  (pl/ Pu  =  2) 
are  left  to  evolve  between  two  planes  (y  =  -1  m  and  y  =  +1  m),  with  gravity  oriented 
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in  the  upward  direction  (g  =  (0,+l}T).  The  boundaries  are  adiabatic  solid  walls. 
The  remainder  of  the  fluid  initial  conditions  above  and  below  the  diaphragm,  U  &  L 
respectively,  are  found  in  Table  4.1.  The  solutions  produced  by  the  MP5,  AW5,  and 
MW5  scheme  were  compared  for  this  problem  against  the  test  solution  at  various 
times  in  Figures  4. 6-4. 9.  From  these  figures,  it  becomes  quite  clear  that  sharp  features 
and  rolled  up  vortices  in  the  MW5  solution  are  far  superior  to  that  of  the  solutions 
produced  by  MP5  and  AW5. 


4.2  Ideal  Magnetohydrodynamics(MHD) 

The  ideal  MHD  test  cases  will  demonstrate  that  the  physical  waves  (entropy,  fast  & 
slow  magneto-acoustic,  and  Alfven)  are  properly  captured  under  various  configura¬ 
tions  by  the  present  schemes. 


4.2.1  ID  MHD  Shock  Tube  Problems 


The  Brio-Wu  problem  [60]  was  used  to  ensure  the  numerical  scheme  sufficiently 
captured  all  of  the  important  features,  ie,  a  contact  discontinuity,  fast  shock,  fast 
rarefaction,  compound  wave,  and  slow  shock.  The  initial  conditions  are  analogous 
to  Sod’s  shock  tube  problem  where  the  initial  left,  and  the  right  states  are  as  follows: 


(pZj  U'Xjli  U y,li  V'Zjh  Bx,li  By^l ,  Bz  iy  Pi ) 

(pr,  Ux,n  ^y,ri  ^ z,ri  Bx,ri  By,r,  Bzr}  Pr) 


=  (0.1,  0,0,  0,0.75, -1,0,1)  x  <  0 
=  (1,0,  0,0, 0.75,  +1,0, 10)  £>0 


(4.4) 


where  7  =  2  and  x  G  (—1, 1].  Figure  4.10  -  4.15  show  the  distributions  of  various 
fluid  properties  at  the  solution  time,  t  =  0.1.  Figure  4.11(a)  shows  a  zoom  in 
of  the  compound  wave  at  x  ~  —.03m.  From  this  figure,  it  is  clear  that  the  MW5 
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scheme  is  able  to  resolve  the  major  features  (compound  wave,  fast  &  slow  rarefaction) 
better  than  MP5  &  AW5.  All  of  the  schemes  (MP5,  AW5,  &  MW5)  resolve  the 
contact  discontinuity  fairly  well,  shown  in  Figure  4. 11(b). The  undershoot  observed 
in  the  velocity  profile,  Figure  4.12(a),  at  x  ^  0.35  m  is  a  well  documented  feature 
in  literature  for  higher  order  MHD  schemes  at  magneto-sonic  points;  Jiang  and  Shu 
[61]  suggest  performing  the  test  problem  in  a  moving  reference  frame  in  order  to 
suppress  the  oscillations. 

4.2.2  Orszag-Tang  Problem 

Various  2D  MHD  test  cases  were  performed  to  ensure  that  the  solver  correctly  cap¬ 
tured  all  physical  waves  and  to  ensure  the  solution  remains  divergence  free,  VB  =  0. 
The  divergence  cleaning  procedure  used  when  performing  MHD  simulations  is  given 
in  Appendix  C.  The  first  numerical  test  case  is  that  of  the  Orszag-Tang  vortex  prob¬ 
lem,  first  introduced  by  Orszag  &  Tang[62],  This  is  a  well-known  model  problem  for 
testing  the  transition  to  supersonic  2D  MHD  turbulence.  The  initial  conditions  of 
the  problem  are  given  by: 


(p,  ux,  uy,  Bx,  By ,  P)  =  (2.778,  -sin (y),  sin(x),  -sin(y)y7pj,  sin(2x)v//pj,  1.667) 

(4.5) 

where  7  =  |.  The  problem  is  set  on  a  periodic  domain  with  the  dimensions  x  :  y  G 
[0,  2n)m  :  [0,  2n)m.  Figure  4.16  shows  the  temperature  distribution  of  the  solutions 
at  t  =  3s.  Although  the  images  illustrated  appear  similar,  upon  further  inspection 
the  flow  features  in  Figures  4.16(a)  &  (c)  are  sharper  than  those  in  Figure  4.16(b). 
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4.2.3  Balsara’s  Rotor  Problem 

The  following  test,  which  was  original  introduced  by  Balsara  &  Spicer [63],  consists  of 
the  propagation  of  strong  torsional  Alfven  waves  into  an  ambient  fluid.  The  problem 
consists  of  a  dense,  rapidly  spinning  cylinder  of  fluid(the  rotor)  surrounded  by  a 
light,  stationary  fluid  (ambient  fluid).  Because  there  is  no  diffusive  transport  in  this 
problem,  the  two  fluids  are  connected  by  an  initially  uniform  magnetic  field.  The 
rapidly  spinning  rotor  causes  torsional  Alfven  waves  to  propagate  into  the  ambient 
fluid,  which  will  lead  to  a  decrease  of  angular  momentum  in  the  rotor.  The  magnetic 
field  is  strong  enough  that  as  it  wraps  itself  around  the  rotor,  the  increased  magnetic 
pressure  will  compress  the  rotor  into  an  oblong  shape.  Balsara  &  Spicer  applied  a 
slight  taper  to  the  initial  density  and  velocity  of  the  rotor  as  to  avoid  generating 
strong  start-up  transient  from  the  computational  scheme.  The  computation  domain 
is  described  by  x  :  y  G  [0, 1]  :  [0, 1]  and  the  initial  conditions  are  described  in  Table 
4.2. 

At  the  solution  time,  t  =  0.295  s,  MW5  and  MP5  shown  in  Figures  4.17  &  4.18, 
respectively,  agree  quite  well  with  the  results  of  Balsara  &  Spicer [63]  using  CFL  = 
0.3,  while  AW5  is  unstable  at  this  CFL  number.  Figure  4.19  illustrates  that  AW5 
can  remain  stable  and  resolve  Balsara’s  rotor  problem  quite  well  by  reducing  the 
CFL  to  0.15. 

4.2.4  Rayleigh  Taylor  MHD  Instability 

Equally  as  important  as  the  stability  of  the  chemical  processes  is  that  of  the  stability 
of  the  MHD.  As  previously  demonstrated  in  Section  4.1.5,  Rayleigh- Taylor  Instabil¬ 
ities  (RTI)  can  arise  from  infinitesimal  disturbances  in  amplitude  and  grow  because 

71 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


of  the  gradient  in  the  applied  force,  e.g.,  a  buoyant  force.  In  [59],  Chandrasekhar  an¬ 
alytically  demonstrated  RTI  with  a  uniform  tangential  magnetic  held  in  both  fluids, 
where  the  stability  growth  rate  is: 


n 


2 


gk 


P 2  -  Pi 
P2  +  Pi 


B 2k2cos2 9 
2n(p2  +  pi) 


(4.6) 


where  k  is  the  wave  number,  n  is  the  growth  rate,  and  9  is  the  angle  between  B 
and  k.  The  critical  strength  of  a  magnetic  held  to  suppress  instability  of  a  mode  of 
wavelength  A  is: 


B„  = 


V9KP2  -pi) 


cos  9 


(4.7) 


Similarly,  the  critical  wavelength  for  a  given  strength  of  a  magnetic  held  can  be 
expressed  as: 


B2cos2  9 
g(p 2  -  Pi) 


(4.8) 


where  A  <  Ac  are  suppressed. 


As  specihed  by  Remade  et  al.  [64],  a  domain  with  dimensions  x  :  y  G  [0,  0.25]m  : 
[— 0.5,0. 5]m  is  enclosed  by  rehective,  adiabatic  walls.  The  heavy,  upper  fluid  is 
separated  from  the  light,  lower  fluid  at  y  —  0.01cos(8-7nr).  The  acceleration  due  to 
gravity  is  g  =  {0,  —  1}T.  The  initial  conditions  conditions  of  the  fluid  are  listed  in 
Table  4.1.  Figures  4.20(a)  -  (c)  show  the  solution  of  the  density  distribution  of  MW5, 
MP5,  and  AW5,  respectively,  at  t  =  2  s  with  grid  resolution  Ax  =  Ay  =  Due 

to  the  lack  of  rolled  up  vortices  in  Figure  4.20(c),  AW5  cannot  properly  resolve  the 
instability  of  the  current  problem.  Thus,  only  MW5  and  MP5  will  be  used  in  the 
next  portion  of  the  test. 


The  suppression  of  the  instability  growth  with  tangential  magnetic  held  strengths 
of  Bx  =  0,  0.2 Bc,  0.5 Bc,  and  0.8 Bc  at  t  =  2  s  using  MW5  &  MP5,  respectively,  is 
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illustrated  in  Figures  4.21  &  4.22  .  The  illustrations  support  the  well  established 
theory  of  Chandrasekhar [5 9]  that  there  is  a  critical  tangential  magnetic  field  that 
will  suppress  a  perturbations  growth  in  this  particular  test  case.  A  normal  magnetic 
field  with  the  strengths  By  =  0.5i?c,  1.5/?c,  and  2 Bc  are  illustrated  at  t  =  2  s  in 
Figures  4.23  &  4.24.  The  normal  magnetic  field  did  not  have  a  significant  affect  on 
the  growth  of  the  instabilities  until  a  larger  field  was  applied  relative  to  the  tangential 
field  strength,  which  is  confirmed  by  Jun  et  al.  [65] .  Additionally,  Jun  et  al.  show 
that  the  growth  of  the  instability  in  the  nonlinear  regime  is  enhanced  by  the  normal 
field  up  to  a  certain  field  strength. 

4.3  Two-Temperature(2T)  Model 

When  there  is  significant  ionization  due  to  shocks  or  non-isentropic  processes,  the 
electron  temperature,  Te,  is  adiabatic  and  must  relax  to  the  translational  temperature 
of  the  heavy  particles,  T^.  Next,  we  present  a  ID  test  case  of  a  Mach  10  fully 
ionized  plasma  of  argon  passing  through  a  normal  shock.  Figure  4.25(a)  illustrates 
the  electron  and  heavy  particle  temperatures  when  they  are  conserved  separately 
and  are  not  allowed  to  relax  toward  equilibrium,  while  Figure  4.25(b)  illustrates  the 
electron  and  heavy  particle  temperatures  relaxing  via  electron-ion  collisions  toward 
equilibrium.  In  these  figures  we  see  that  MP5  is  able  properly  conserve  electron 
entropy  as  well  relax  to  the  correct  equilibrium  temperature,  Te  =  Th  =  7800/1 . 
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upper  part 

lower  part 

p 

2 

1 

p 

2-2  y 

2 -y 

'LL X 

exsin(8'7rx)cos(7ry)sinT~1(7ry ) 

same  as  upper  part 

Uy 

— eycos(8'Kx)sinT  (ny) 

same  as  upper  part 

Table  4.1:  RTI  problem  hydrodynamic  initial  conditions  as  specified  by  Remade  et 
al. [64]  where  r  =  6,  M0  =  0.1,  ey  =  M0^/ y/2,  and  ex  =  —  eyr/16. 


r  <  r0 

r0  <  r  <  ri 

r  >  r  1 

P 

10 

1 

1  +  9/ 

P 

0.5 

0.5 

0.5 

vo(y-yo) 

ro 

x  vo(y—yo) 

J  r0 

0 

Uy 

v0  ( x—xo ) 
ro 

rv  0(x—xo) 

J  r0 

0 

Table  4.2:  Balsara’s  rotor  problem  initial  conditions  as  specified  by  Toth[66]  with  a 
magnetic  field  Bx  =  2.5//i0,  where  r0  =  0.1,  rx  =  0.115,  /  =  v0  =,  n0  =  1, 

(xo,yo)  =  (0.5, 0.5),  and  7  =  5/3. 
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(b) 

Figure  4.1:  Density  distribution  of  Sod’s[54]  ID  Shock  Tube  problem  at  runtime  = 
2  s  and  CFL=0.4,  where  (a)  compares  the  MP5,  AW5,  and  MW5  solutions  with 
Ax  =  100/im  and  (b)  is  an  expanded  view  near  the  contact  discontinuity  showing 
the  convergence  of  MP5  with  successively  reduced  grid  size  Ax. 
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Figure  4.2:  Density  distribution  of  Woodward  and  Collela’s[55]  ID  Blastwave  prob¬ 
lem  using  MP5,  AW5,  and  MW5,  Ax  =  2.5  x  10”3m,  CFL  =  0.4,  runtime  =  38ms. 
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Figure  4.3:  Density  distribution  of  the  contact  discontinuities  for  the  ID  Blastwave 
problem  [55]  using  MP5,  AW5,  and  MW5,  Ax  =  2.5  x  10_3m,  CFL  =  0.4,  runtime 
=  38ms,  where  (a)  and  (b)  are  expanded  views  near  contact  discontinuities  of  Figure 
4.2. 
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5 


Figure  4.4:  Density  distribution  of  the  ID  Shock-Entropy  Interaction  problem  [56] 
using  MP5,  AW5,  and  MW5,  Ax  =  6.67  x  10_3m,  CFL  =  0.4,  runtime  =  360ms, 
where  (a)  is  the  wide  view  and  (b)  is  an  expanded  view  near  an  entropy  distrubance. 
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Figure  4.5:  Mach  2.4  flow  over  a  backward  facing  step  solution  (left)  using  MW5 
with  Ax  =  Ay  =  pp4m  compared  to  the  experimental  results  (right)  [57]. 
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(a)  MP5  (b)  AW5  (c)  MW5  (d)  High  Res. 

Figure  4.6:  Density  distribution  for  hydrodynamic  Rayleigh  Taylor  Instability  prob¬ 
lem  at  t  —  0.75s,  where  Ax  =  Ay  =  ^ m  (For  High  Resolution,  MW5  was  used 
with  Ax  =  Ay  =  6.25  x  10_4m),  min  ■  max. 
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(a)  MP5  (b)  AW5  (c)  MW5 


(d)  High  Res. 


Figure  4.7:  Density  distribution  for  hydrodynamic  Rayleigh  Taylor  Instability  prob¬ 
lem  at  t  —  1.50s,  where  Ax  =  Ay  =  ^ m  (For  High  Resolution,  MW5  was  used 
with  Ax  =  Ay  =  6.25  x  10_4m),  min  ■  max. 
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(a)  MP5  (b)  AW5  (c)  MW5  (d)  High  Res. 

Figure  4.8:  Density  distribution  for  hydrodynamic  Rayleigh  Taylor  Instability  prob¬ 
lem  at  t  —  2.25 s,  where  Ax  =  Ay  =  ^ m  (For  High  Resolution,  MW5  was  used 
with  Ax  =  Ay  =  6.25  x  10_4m),  min  ■  max. 
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Figure  4.9:  Density  distribution  for  hydrodynamic  Rayleigh  Taylor  Instability  prob¬ 
lem  at  t  —  3.00s,  where  Ax  =  Ay  =  ^ m  (For  High  Resolution,  MW5  was  used 
with  Ax  =  Ay  =  6.25  x  10_4m),  min  ■  max. 
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Figure  4.10:  Density  distribution  of  the  Brio-Wu[60]  problem  using  MP5,  AW5  and 
MW5,  Ax  =  2.5  x  10~3m,  CFL  =  0.4,  runtime  =  100ms 
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Figure  4.11:  Density  distribution  for  the  ID  Brio-Wu[60]  problem  using  MP5,  AW5, 
and  MW5,  Ax  =  2.5  x  lCT3m,  CFL  =  0.4,  runtime  =  100ms,  which  illustrate  the 
expanded  views  near  the  compound  wave  (a)  and  contact  discontinuity  (b)  of  Figure 
4.10. 
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Figure  4.12:  Velocity  distribution  of  the  x-component  for  the  ID  Brio-Wu[60]  prob¬ 
lem  using  MP5,  AW5,  and  MW5,  Ax  =  2.5  x  10_3m,  CFL  =  0.4,  runtime  =  100ms, 
which  (b)  illustrates  the  expanded  views  near  the  trailing  edge  of  the  fast  rarefaction 
of  (a). 
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Figure  4.13:  Velocity  distribution  of  the  y-component  for  the  Brio-Wu[60]  problem 
using  MP5,  AW5  and  MW5,  Ax  =  2.5  x  10_3m,  CFL  =  0.4,  runtime  =  100ms 
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Figure  4.14: 
MW5,  Ax  = 


Pressure  distribution  of  the  Brio-Wu[60]  problem  using  MP5,  AW5  and 
2.5  x  lCT3m,  CFL  =  0.4,  runtime  =  100ms 
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Figure  4.15:  By  distribution  of  the  Brio-Wu[60]  problem  using  MP5,  AW5  and  MW5, 
Ax  =  2.5  x  lCT3m,  CFL  =  0.4,  runtime  =  100ms 
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Figure  4.16: 
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Figure  4.17:  The  density,  thermal  pressure,  Mach  number,  and  the  magnetic  pressure 
distributions  at  t  =  0.295s  for  Balsara’s  rotor  problem.  The  solution  was  obtained 
using  MW5  with  a  grid  resolution  of  Ax  =  Ay  =  -j and  CFL  =  0.3.  The  30 
contour  lines  are  shown  for  the  ranges  0.532  <  p  <  10.83^|,  .007  <  P  <  0.702  Pa, 
0  <  M  <  3_g4;  and  0.007  <  <  0.702 Pa,  as  prescribed  by  Toth  [66]. 
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Figure  4.18:  The  density,  thermal  pressure,  Mach  number,  and  the  magnetic  pressure 
distributions  at  t  =  0.295s  for  Balsara’s  rotor  problem.  The  solution  was  obtained 
using  MP5  with  a  grid  resolution  of  Ax  =  Ay  =  and  CFL  =  0.3.  The  30 

contour  lines  are  shown  for  the  ranges  0.532  <  p  <  10.83^|,  .007  <  P  <  0.702  Pa, 
0  <  M  <  3_g4;  and  0.007  <  <  0.702 Pa,  as  prescribed  by  Toth  [66]. 
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Figure  4.19:  The  density,  thermal  pressure,  Mach  number,  and  the  magnetic  pressure 
distributions  at  t  =  0.295s  for  Balsara’s  rotor  problem.  The  solution  was  obtained 
using  AW5  with  a  grid  resolution  of  Ax  =  Ay  =  and  CFL  =  0.15.  The  30 
contour  lines  are  shown  for  the  ranges  0.532  <  p  <  10.83^|,  .007  <  P  <  0.702  Pa, 
0  <  M  <  3_g4;  and  0.007  <  <  0.702 Pa,  as  prescribed  by  Toth  [66]. 
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Figure  4.20:  Growth  of  instability  without  the  presence  of  a  magnetic  field  at  t  =  2s 
for  MW5,  MP5,  and  AW5  where  Ax  =  Ay  =  Density  distribution  with  20 

density  contours  between  1  |||  2. 
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(a)  Bx  =  0.0 Bc 


Br  =  0.2  Br 


Figure  4.21:  The  impact  of  the  tangential  magnetic  field,  Bx,  on  the  growth  of 
instabilities  at  t  —  2s  where  Ax  =  Ay  =  ~^m  and  Bc  =  0.5T  using  MW5.  Density 
distribution  with  20  density  contours  between  1  1 1  ■  2. 
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Figure  4.22:  The  impact  of  the  tangential  magnetic  field,  Bx,  on  the  growth  of 
instabilities  at  t  =  2s  where  Ax  =  Ay  =  and  Bc  =  0.5T  using  MP5.  Density 
distribution  with  20  density  contours  between  1  1 1  |  2. 
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(a)  By  =  0.5 Bc  (b)  By  =  l.0Bc  (c)  By  =  1.5 Bc  (d)  By  =  2.0 Bc 


Figure  4.23:  The  impact  of  the  normal  magnetic  field,  Bx,  on  the  growth  of  insta¬ 
bilities  at  t  =  2s  where  Ax  =  Ay  =  and  Bc  =  0.5T  using  MW5.  Density 
distribution  with  30  density  contours  between  1  1 1  |  2. 
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(a)  By  =  0.5 Bc  (b)  By  =  l.0Bc  (c)  By  =  1.5 Bc  (d)  By  =  2.0 Bc 


Figure  4.24:  The  impact  of  the  normal  magnetic  field,  Bx,  on  the  growth  of  insta¬ 
bilities  at  t  =  2s  where  Ax  =  Ay  =  and  Bc  =  0.5T  using  MP5.  Density 

distribution  with  30  density  contours  between  1  m  2. 
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Figure  4.25:  Mach  10  fully  ionized  Argon  shock  in  shock  reference  frame  with(a)  and 
without (b)  Te  relaxation  where  Ax  =  using  MP5. 
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CHAPTER  5 


Simplified  Approach  for  PDRIME  Simulations 


Now  that  the  various  numerical  schemes  have  been  validated,  we  first  investigate 
the  performance  of  a  PDRE  flight  configuration,  and  later  the  PDRIME,  using  a 
simplified  non-reactive  model,  Equation  2.5.  In  the  simplification,  we  assumed  that 
all  reactants  ( H2  and  O2 )  in  the  combustion  chamber  have  been  consumed,  leaving 
just  products  ( H20 )  and  the  seeded  Cesium.  We  also  assumed  that  the  detonation 
has  left  the  chamber.  This  model  has  been  shown  to  replicate  blowdown  conditions 
in  a  rocket  nozzle  reasonably  accurately,  as  shown  in  Figure  5.1  from  the  work  by 
Cambier.  Following  Cambier’s  blowdown  model[67],  we  initialize  the  combustion 
chamber  with  the  post-combustion  and  post-detonation  conditions: 


Po  =  100  atm 
f0  =  3000  K 


P  0  = 


-  _Ql 


RT0 


Rather  than  simulating  the  combustion  chamber,  the  combustion  chamber  is  mod¬ 
eled  as  a  “reservoir”  whose  conditions  (Po,  To)  temporally  evolve  as  prescribed  by 
Cambier’s  blowdown  equations[67]: 


P0  =  P0[/(t)]7/(7-1}  (5.1a) 
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po  =Po[f(t)}1/b  1} 


(5.1b) 


T0  =  T0[f(t)]  (5.1c) 

where  f(t)  =  and  v  the  blowdown  frequency  and  is  equal  to  87.5  Hz  for 
this  particular  case.  Now,  by  using  the  isentropic  relations,  we  can  determine  the 
conditions  at  the  PDRE  throat  as  a  function  of  time.  With  this  simplification,  our 
computational  domain  no  longer  needs  to  include  the  combustion  chamber,  this  is 
illustrated  in  Figure  5.2.  Next,  we  assume  that  the  back  pressure,  the  atmospheric 
pressure  Patm,  is  sufficiently  low,  so  that  there  is  no  back  flow  into  the  nozzle.  Pres¬ 
sures  at  the  altitudes  we  are  simulating  are  high  enough  that  even  with  the  expan¬ 
sion  in  the  nozzle,  the  nozzle  exit  pressure  is  greater  than  the  atmospheric  pressure, 
Pexit  >  Patm ■  Because  of  the  extremely  high  flight  altitudes  (altitude  ~  15  —  20 
km)  ,  the  PDRE  performance  does  not  have  any  significant  variance  with  regard  to 
altitude. 

Next,  we  extend  this  simplified  model  of  the  PDRE  to  the  PDRIME  with  a 
‘Magnetic  Piston’.  In  this  extension,  we  assume  that  if  the  temperature  is  greater 
than  3000K;  the  conductivity  is  a  constant  non-zero  value,  in  this  particular  case 
er  =  1000 mho/m  in  the  nozzle  and  a  =  500 mho/m  in  the  bypass  tube.  The  nozzle 
impulse  contribution  is  calculated  the  same  way  as  the  PDRE,  while  bypass  tube 
impulse  contribution  is  calculated  by  integrating  the  Lorentz  force  over  the  cycle 
duration,  as  follows: 

hypass  =  [  (j  x  B)x  dr  «  f  auBz(Kx  -  1)  dr  (5.2) 

Jo  Jo 
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where  all  variables  have  been  previously  defined.  We  also  assumed  a  constant  mag¬ 
netic  held  for  both  the  nozzle  and  bypass  section  of  the  PDRIME,  as  illustrated  in 
Figure  1.7.  The  configuration  called  ‘NG’  refers  to  the  PDRE  with  energy  extraction 
in  the  nozzle  (for  use  of  the  ‘nozzle  generator’),  so,  as  expected  there  is  a  drop  in 
impulse  as  shown  in  Figure  5.3.  A  promising  concept  postulated  by  Cambier[40] 
involving  the  PDRE  with  ‘Magnetic  Piston’  was  also  simulated,  but  due  to  the  lack 
of  energy  supplied  to  the  ‘Magnetic  Piston’  by  the  generator,  this  flight  configuration 
could  not  be  effectively  utilized.  In  Chapter  7,  the  true  energy  ‘cost’  of  the  chamber 
piston  will  be  investigated  by  exploring  the  nature  of  the  ionized  gas  and  chemical 
kinetics  subjected  to  the  Lorentz  force.  We  will  now  explore  other  PDRE  flight  con¬ 
figurations  which  will  build  off  of  the  PDRE  with  ‘NG’  concept  in  order  to  further 
optimize  the  PDRE. 

5.1  PDRIME  with  Bypass  Configurations 

There  were  multiple  flight  configurations  as  well  as  different  flight  condition  tested. 
Using  the  standard  combustion  chamber  condition,  nozzle,  and  nozzle  generation 
configuration  we  were  able  to  test  various  configurations  of  the  Bypass  Tube.  Figure 
5.3  shows  the  performance  of  the  standard  PDRIME  configuration  under  various 
flight  condition  as  they  compare  the  performance  of  a  PDRE,  which  we  will  refer  to 
as  the  baseline.  With  this  configuration,  there  are  some  marginal  performance  gains 
above  the  baseline  case.  The  standard  PDRIME  utilizes  a  constant  magnetic  held, 
which  is  active  only  if  the  fluid  temperature  is  greater  than  3000K  and  the  fluid  is 
moving  in  the  stream-wise  direction.  When  activated,  a  constant  magnetic  held  of 
3T  is  applied  to  the  applicable  fluid. 
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5.2  Temporal/Temperature  Controllers  of  the  Magnetic  Field 


We  performed  component  evaluations  for  the  various  flight  configurations  with  dif¬ 
ferent  assumptions.  This  involved  examining  the  net  energy  consumed  compared  to 
impulse  gained  or  lost  by  a  particular  component  of  the  flight  configuration.  During 
the  course  of  component  performance  evaluation,  it  was  determined  that  an  MHD 
bypass  accelerator,  where  a  “piston”  like  acceleration  was  applied  in  the  bypass  sec¬ 
tion  as  done  for  the  magnetic  chamber  piston,  was  under-performing  relative  to  the 
nozzle  which  is  shown  in  Figure  5.4.  Zeineh[33]  also  saw  the  same  under-performance 
in  his  2-D  simulations.  The  bypass  available  energy,  A Egen,  is  being  applied  to  the 
bypass  over  approx  0.1  ms  within  a  relatively  small  volume,  dV,  wherein  most  of  the 
available  energy  is  directed  towards  Joule  heating  rather  than  accelerating  the  fluid. 
This  leads  to  large  temperature  spikes  and  little  impulse  increase. 

From  E  =  fcv  dT  +  \pv2  with  introduction  of  Energy  at  constant  volume 

AEgen  =  cvAT  +  T^pAv2  (5.3) 

We  can  see  where  and  how  the  energy  is  being  used  in  the  bypass.  AEgen  is  a  fixed 
amount  of  energy  produced  by  the  MHD  generator  in  the  nozzle.  cvAT  represents  the 
amount  of  energy  that  is  converted  to  internal  energy,  which  for  our  simplified  model 
of  constant  heat  capacity  is  represented  by  a  rise  in  temperature.  Lastly,  \pAv2  oc 
AKE  represents  the  amount  of  kinetic  energy  which  is  imparted  upon  the  fluid,  and 
since  energy  is  introduced  at  constant  volume,  p  remains  constant,  so  this  change  in 
kinetic  energy  is  directly  related  to  our  change  in  impulse,  I  oc  AKE^.  In  order  to 
reduce  the  temperature  spikes  and  thus  increase  the  impulse,  when  applying  MHD  to 
the  bypass,  one  must  regulate  the  AEgen  available  for  use  by  the  bypass  accelerator. 

To  control  this  energy  reintroduction,  one  must  either  control  the  conductivity  of  the 
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fluid,  a,  or  the  strength  of  the  magnetic  held,  B ,  acting  on  the  bypass.  Since  it  is 
impractical  if  not  impossible  to  dynamically  control  the  conductivity  in  the  bypass, 
we  chose  the  latter. 

Our  first  approach  is  to  scale  the  maximum  magnetic  held  in  the  bypass  with 
time,  for  which  we  prescribe  a  simple  algorithm.  First,  the  magnetic  held  is  initially 
zero  until  the  initiation  time,  to,  when  at  least  one  grid  cell  in  the  bypass  meets  the 
normal  criteria  for  MHD  accelerator:  T,  >  3000/i  and  ?i,  >  0.  The  magnetic  held  is 
then  scaled  with  the  maximum  held  applied,  B0.  The  magnetic  held  in  the  bypass 
is  prescribed  as  follows: 


(  0 


B(t )  —  n(t)  •  Bq  to  <  t  <C  to  T  At 
^  B0  t  >  to  +  At  J 


(5.4) 


B0 :  Maximum  Magnetic  Field  Strength 
t0:  time  at  which  T  >  3000 K  and  u  >  0 
At:  ramp  up  time 
fl(f):  Magnetic  Field  Scaling  Factor 


We  ran  a  series  of  simulations  incorporating  this  prescribed  magnetic  held,  adjust¬ 
ing  both  fl(f)  and  At  as  a  means  of  “open  loop”  control.  Our  tests  utilized  three 
different  magnetic  held  scaling  factors: 

11(f)  =  linear  scaling  where  B  oc  t 


11(f)  =  2  energy  scaling  ( Energy  oc  B 2),  where  B  oc  O 

II (t)  =  f)n  power  law  scaling  (n  =  3)  where  B  (xtn 
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Figure  5.5  compares  the  evolution  of  these  magnetic  held  scaling  factors  with  time. 
Figure  5.6  shows  that  the  temporally  varied  bypass  magnetic  held  with  the  “open 
loop”  control  shown  in  Equation  5.4  has  the  same  marginal  performance  improve¬ 
ments  performance  as  the  standard  bypass  configuration.  The  optimal  configuration 
of  the  bypass  magnetic  held  is  dictated  by  the  flight  conditions  which  include  flight 
Mach  number  and  altitude. 

In  our  next  approach,  we  return  to  the  same  problem  of  controlling  the  rate 
at  which  energy  is  consumed  by  the  bypass.  It  was  previously  stated,  that  it  was 
impractical  to  dynamically  control  the  conductivity  in  the  bypass,  a.  When  evalu¬ 
ating  the  nature  of  the  conductivity  of  a  fluid  in  the  current  regime,  Equation  2.30 
shows  the  conductivity’s  strong  dependency  on  the  temperature,  T.  Rather  than 
dynamically  controlling  the  conductivity  of  the  fluid,  we  modeled  the  behavior  of  an 
‘ideally’  conductive  fluid  by  prescribing  its  dependency  to  temperature,  a  =  cr(T). 
We  wish  to  implement  this  “closed  loop”  controller  in  the  bypass  in  order  to  scale 
the  magnetic  held,  Bi,  with  the  sensible  local  temperature  in  the  bypass,  T%.  The 
normal  criteria  for  MHD  acceleration  are  still  utilized.  The  bypass  magnetic  held 
“feedback”  function  is  prescribed  as  follows: 


(  0  T  <  Tmin 

B(T )  =  0(T)  •  B0  Tmin  <T<  (Tmax  -  Tref) 

^  Bq  T  >  (Tmax  —  Tref) 

Bq:  Maximum  Magnetic  Field 
Tmin:  Minimum  Temperature  (3000/v) 

0(T):  Magnetic  Field  Scaling  Factor 
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@(T)  =  ( T 1  ^  { )  linear  scaling  B  ocT 
Tref  and  Tmax  are  adjusted  to  optimize  performance 

The  term  “feedback”  here  is  used  in  the  sense  that  temperature  is  measured,  and  the 
magnetic  field  is  adjusted  according  to  Equation  5.5.  Figures  5.7  through  5.10  show 
that  the  “closed-loop”  controller  in  the  bypass  tube  is  able  to  match  the  performance 
of  the  standard  configuration  at  a  flight  altitude  of  25km  and  Tmax  <  8  x  103Ji . 
But  at  other  flight  altitudes,  this  “closed  loop”  controller  falls  short.  The  B{T ) 
model  is  much  more  feasible  for  implementing  into  a  physical  system  because  of  the 
temperature  dependency  of  many  of  the  parameters  that  contribute  to  impulse  from 
the  bypass  accelerator,  e.g.,  ct(T)  and  B{T). 

5.3  PDRIME  with  2D  Bypass  Configuration 

There  was  further  examination  of  the  PDRIME  with  various  configurations  per¬ 
formed  in  Zeineh[33]  and  Zeineh  et  al. [34].  In  their  study,  Zeineh  et  al.  found  the 
flight  configuration  of  the  PDRIME  at  altitudes  20,  25,  and  30  km  at  relatively 
low  Mach  numbers  (M  <  5)  had  significant  performance  increases  over  the  base¬ 
line  PDRE.  Figures  5.11  through  5.13  show  some  of  the  performance  gains  of  the 
PDRIME  at  different  altitudes  with  the  optimal  bypass  section  lengths(Lfeypass  =  3, 
4,  and  6  m)  demonstrated  in  [34],  While  the  results  of  Zeineh’s  multidimensional 
simulations  were  quite  promising,  the  effects  of  complex  kinetics  (i.e.  hydrogen-air 
chemistry  and  ionization  process)  were  not  rigorously  investigated.  The  constant 
conductivity  assumption  does  not  allow  for  the  “closed  loop”  evolution  of  the  con¬ 
ductivity  (i.e.  Joule  heating  further  ionizing  fluid),  therefore  in  Chapter  7  we  will 
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discuss  and  examine  the  effects  of  MHD  when  the  conductivity  is  described  more 
accurately  via  Equation  2.27. 
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Time  [ms] 


Figure  5.1:  Variation  in  impulse  for  a  PDRE.  Results  are  shown  from  a  full  quasi-ID 
transient  PDRE  simulation  and  a  cycle  approximated  by  a  constant  volume  reaction 
and  a  blow-down  period(from  Cambier[67]) 
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, Bypass  inlet  B.C. 


7 


Nozzle  throat  B.C. 

Figure  5.2:  Quasi  ID  Computational  Domain,  where  the  Bypass  and  Nozzle  inlet 
and  exit  boundary  conditions  are  shown 
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Figure  5.3:  Quasi  ID  Performance:  Impulse  Loss  in  Nozzle  Energy  Generation  with 
H20  product 
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Figure  5.4:  Quasi  ID  Component  Performance:  Bypass  vs.  Nozzle  Impulse  at  M  = 
9  and  Alt:  25km  with  H20  product 
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Figure  5.5:  Magnetic  Field  Strength,  B,  as  a  function  of  time,  using  various  Magnetic 
Field  Strength  functions,  11(f) 
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Figure  5.6:  PDRIME:  effects  of  flight  Mach  number.  Magnetic  Field,  B(t)  ~  t2 
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Figure  5.7:  PDRIME  with  Chamber  Piston:  effects  of  flight  rnach  number.  Magnetic 
Field  B(T),  Tre f  =  OK  and  Tmax  =  6  x  10 3K 
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Figure  5.8:  PDRIME:  effects  of  flight  mach  number.  Magnetic  Field  B(T),  Tre{  =  OK 
and  Tmax  =  7  x  103A" 


115 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


Figure  5.9:  PDRIME:  effects  of  flight  mach  number.  Magnetic  Field  B(T),  Tre{  =  OK 
and  Tmax  =  8  x  103A" 
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Figure  5.10:  PDRIME:  effects  of  flight  mach  number.  Magnetic  Field  B(T), 
Tre f  =  0A"  and  Tmax  =  9  x  103A" 
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Figure  5.11:  PDRIME  impulse  per  cycle  at  20  km  plot  against  various  Mach  numbers 
and  bypass  area  per  unit  depth.  The  chamber  is  initially  seeded  with  0.5%  cesium 
by  number  at  an  initial  temperature  of  3000  K.  The  bypass  length  is  Lbypass  =  3  nr 
and  is  seeded  with  0.1%  cesium  by  number,  (from  Zeineh  et  al.  [34]) 
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Figure  5.12:  PDRIME  impulse  per  cycle  at  25  km  plot  against  various  Mach  numbers 
and  bypass  area  per  unit  depth.  The  chamber  is  initially  seeded  with  0.5%  cesium 
by  number  at  an  initial  temperature  of  3000  K.  The  bypass  length  is  Lbypass  =  4  m 
and  is  seeded  with  0.1%  cesium  by  number,  (from  Zeineh  et  al.  [34]) 
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Figure  5.13:  PDRIME  impulse  per  cycle  at  30  km  plot  against  various  Mach  numbers 
and  bypass  area  per  unit  depth.  The  chamber  is  initially  seeded  with  0.5%  cesium 
by  number  at  an  initial  temperature  of  3000  K.  The  bypass  length  is  Lbypass  =  6  m 
and  is  seeded  with  0.1%  cesium  by  number,  (from  Zeineh  et  al.  [34]) 
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CHAPTER  6 


Detonation  Stability  Phenomena 


This  chapter  was  taken  with  slight  modification  from  the  article  “Stability  of  flame- 
shock  coupling  in  detonation  waves:  ID  dynamics” ,  published  in  the  journal  Com¬ 
bustion  Science  and  Technology  [68]. 

6.1  Ignition  and  Instabilities 

Direct  initiation  of  the  ID  detonation  and  the  ensuing  instabilities  are  now  exam¬ 
ined,  in  part  for  fundamental  understanding  but  also  to  establish  a  baseline  against 
which  the  effects  of  MHD  may  be  compared  in  Chapter  7.  In  this  study,  a  cham¬ 
ber  is  filled  with  a  stoichiometric  mixture  of  H2  and  air  (temperature  300  K  and 
pressure  1  atm),  and  ignition  is  achieved  by  setting  a  region  adjacent  to  an  end- wall 
of  the  simulated  shock  tube  at  high  pressure  (40  atm  for  most  computations)  and 
temperature  (1500  K),  as  a  simulated  spark.  This  direct  initiation  is  preferable  to  a 
deflagration-to  detonation  transition  (DDT),  since  the  latter  is  much  more  sensitive 
to  initial  conditions  and  grid  resolution,  it  requires  inclusion  of  species  diffusion,  and 
it  requires  a  very  long  computational  domain  to  reproduce  both  the  DDT  and  the 
subsequent  evolution  of  the  detonation.  Nevertheless,  even  direct  initiation  is  sen¬ 
sitive  to  initial  conditions  and  resolution.  The  requirements  to  achieve  detonation 
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ignition  with  the  MP5  scheme  include  a  grid  cell  size  Ax  of  less  than  50  /mi  and 
a  distributed  simulated  spark  region  (of  length  £spark  ranging  from  0.25  to  0.5  cm) 
with  sufficiently  high  pressure.  For  example,  Figure  6.1(a)  is  an  x-t  diagram  illus¬ 
trating  the  pressure  contours  of  a  spark-ignited  mixture  with  initial  pressure  20  atm, 
^spark  =  0.25  cm,  and  grid  size  Ax  =  50  /irn,  which  does  not  achieve  detonation.  In 
contrast,  Figure  6.1(b)  illustrates  contours  of  the  same  mixture  and  grid  resolution 
but  with  a  higher  spark  pressure,  50  atm,  and  larger  region  (£spark  =  0.5  cm),  which 
does  achieve  detonation.  The  constraints  on  initial  spark  conditions  can  be  related 
to  the  concepts  of  minimum  energy  and  kernel  size  for  direct  initiation  [69];  similar 
studies  by  Eckett  et  al.  [20]  and  He  &  Karagozian[70]  indicate  there  are  both  pressure 
and  temperature  requirements  for  the  spark.  With  an  understanding  of  the  range  of 
satisfactory  kernel/spark  sizes  and  pressures,  we  fix  the  spark  conditions  to  consis¬ 
tently  achieve  a  rapid  initiation,  in  order  to  remove  the  dependence  of  the  long-term 
dynamics  on  the  initial  conditions.  For  the  remainder  of  the  studies  in  this  chapter, 
the  spark  conditions  are  £spark  =  0.25  cm,  Pspark  =  40  atm,  and  Tspark  =  1500  A". 

The  succesful  detonation  initiation  event  proceeds  in  two  phases.  First,  the  gas  in 
the  spark  region  rapidly  burns  and  increases  the  pressure,  in  a  nearly  constant- volume 
combustion  process.  This  high  pressure  generates  a  strong  shock  which  propagates 
into  the  unburnt  mixture,  which  itself  is  ignited  after  a  time  delay  and  rapidly  burns, 
starting  from  the  region  closest  to  the  spark.  In  a  scenario  described  as  the  SWACER 
mechanism  [71],  the  combustion  wave  is  amplified  as  it  overtakes  the  leading  shock, 
and  the  coalescence  of  the  two  fronts  leads  to  extremely  high  peak  pressures  for  a 
very  short  time.  This  event  is  easily  identified  in  the  trace  of  the  peak  pressure 
versus  time  as  shown,  for  example,  in  Figure  6.2,  and  is  referred  to  hereafter  as  the 
“re-explosion”  event,  the  first  explosion  having  taken  place  within  the  initial  spark 
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region.  Two  different  grid  sizes  are  used  in  Figures  6.2(a)  and  6.2(b),  Ax  =  12.5 
fj, m  and  2.5  //m.  respectively,  using  the  MP5  scheme.  A  more  detailed  examination 
of  the  dynamics  for  these  two  different  grid  sizes  is  shown  in  Figures  6.3  and  6.4,  as 
will  be  discussed  below. 

The  high  pressure  of  this  re-explosion  event  initiates  another  strong  shock,  which 
is  followed,  after  an  an  induction  length  (£,  measured  in  the  reference  frame  of 
the  shock)  by  the  combustion  zone.  This  flame  is  initially  strongly  coupled  to  the 
shock  (&  —$■  0)  and  the  wave  is  strongly  over-driven,  i.e. ,  its  speed  exceeds  that  of 
the  Chapman-Jouget  (CJ)  detonation.  As  the  degree  of  overdrive  decays  and  the 
detonation  approaches  the  CJ  limit,  instabilities  begin  to  appear,  as  shown  in  Figure 
6.2(a)  after  about  35  /zs  and  in  Figure  6.2(b)  after  about  25  /zs. 

These  spark-ignited  detonation  simulations  demonstrate  the  appearance  of  differ¬ 
ent  instability  modes.  For  both  sets  of  results  in  Figure  6.2,  instabilities  appear  when 
the  detonation  becomes  close  to  the  CJ  condition,  starting  with  a  small-amplitude, 
but  high-frequency  mode  -  hereafter  referred  to  as  the  ‘HF’  mode.  For  the  low- 
resolution  case  in  Figure  6.2(a),  the  transition  to  the  high-frequency  (HF)  mode 
occurs  at  a  time  of  the  order  of  30  —  35  /zs;  this  instability  regime  is  shown  in  detail 
in  Figure  6.3(a).  For  the  high-resolution  case,  Figure  6.2(b),  the  transition  occurs 
earlier,  close  to  25  /zs;  this  regime  is  shown  in  detail  in  Figure  6.3(b)  for  roughly 
the  same  time  period  as  in  Figure  6.3(a).  For  both  grid  sizes,  we  observe  that  at 
around  45  /zs,  there  is  a  transition  towards  a  lower  frequency  but  high-amplitude 
mode  -  referred  to  here  as  the  ‘HA’  mode.  This  transition  is  more  gradual  in  the 
high-resolution  case,  with  both  modes  coexisting  during  a  period  of  time  (between 
approximately  44  and  48  /zs),  as  shown  in  more  detail  in  Figure  6.4(b).  In  the 

low-resolution  case  of  Figure  6.4(a),  however,  the  behavior  is  less  gradual  and  more 
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chaotic.  The  contrast  between  these  two  profiles  is  striking;  while  a  periodic  signal 
can  still  be  detected  in  Figure  6.4(a),  the  characteristics  and  frequencies  are  both 
very  different.  Besides  the  smoothness  of  the  temporal  waveforms  of  the  high  am¬ 
plitude  instabilities,  we  also  note  that  a  period-doubling  in  the  high  resolution  case, 
i.e.,  the  HA  signal  has  a  dual  oscillation  (high-low  pressures)  which  is  not  apparent 
at  lower  resolution1.  We  note  also  that  there  are  other  manifest  differences  in  the 
specific  dynamical  features  of  these  instabilities,  e.g.,  in  the  appearance  of  noise  in 
the  waveform  after  the  re-explosion  event  (Figure  6.2(a)),  which  is  eliminated  for  the 
higher  resolution  case  in  Figure  6.2(b).  Thus,  it  is  clear  that  there  can  be  significant 
effects  of  the  grid  resolution  on  the  dynamics  of  the  instabilities.  Furthermore,  the 
appearance  of  sharp  features  in  the  traces  also  suggests  that  special  care  must  be 
exercised  in  avoiding  numerical  procedures  which  can  arbitrarily  sharpen  gradients, 
as  mentioned  earlier. 

There  is  no  obvious  “very  high  frequency”  mode  that  arises  after  the  re-explosion 

event  and  before  the  initiation  of  the  high  frequency  mode,  as  seen  by  Leung  et 

al.[17],  in  the  well-resolved  result  in  Figure  6.2(b).  But  if  one  explores  in  detail 

the  time  regime  after  the  re-explosion  event,  for  grid  sizes  Ax  of  2.5  /jrn  and  even 

smaller,  as  shown  in  Figure  6.5,  one  does  observe  relatively  low  amplitude  and  very 

high  frequency  oscillations,  with  frequencies  and  amplitudes  dependent  on  the  grid 

resolution.  To  examine  the  origin  of  these  low  amplitude  oscillations,  the  variation  in 

induction  length  as  a  function  of  time,  determined  from  peaks  in  the  concentration 

of  H  atoms,  may  be  explored  for  these  different  resolutions.  A  plot  of  induction 

length  as  a  function  of  time  is  shown  in  Figure  6.6(a),  with  an  expanded  view  in 

Figure  6.6(b)  corresponding  to  the  same  time  period  as  shown  in  Figure  6.5.  The 
lrThis  was  verified  for  longer  time  periods  than  shown;  for  clarity  purposes,  the  extent  of  the 
simulation  results  shown  in  the  figures  has  been  truncated. 
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results  in  Figure  6.6  indicate  that  the  amplitude  of  oscillations  in  induction  length 
corresponds  to  the  grid  size  Ax,  and  moreover,  the  frequency  of  oscillation  correponds 
to  the  ratio  of  the  CFL  number  here  (0.4)  to  the  sampling  period  (1  nanosecond). 
Thus  it  appears  that  the  very  high  frequency  (and  very  low  amplitude)  oscillations 
seen  using  complex  kinetics  merely  correspond  to  numerical  uncertainties  associated 
with  the  location  of  the  peak  pressure  or  peak  in  atomic  hydrogen  (mole  fraction). 

Because  the  re-explosion  event  is  clearly  identifiable,  we  can  use  this  feature  to 
conduct  a  more  detailed  study  of  the  effect  of  grid  resolution.  For  example,  one  can 
examine  the  variation  of  the  measured  time  delay  to  this  second  explosion  event, 
teXpi  for  the  specific  initial  spark  conditions  noted  previously.  Here  the  uniform  grid 
spacing  Ax  is  varied  from  0.5  /im  to  20  fim.  Figure  6.7  illustrates  how  the  time  to 
re-explosion  for  the  MP5,  AW5,  and  MW5  schemes  varies  with  the  grid  resolution. 
Since  the  AW5  &  MW5  schemes  are  more  diffusive  than  MP5,  it  is  not  surprising 
that  the  AW5  &  MW5  curves  exhibit  a  shallower  profile.  The  most  striking  feature 
here  is  the  non-monotonic  behavior,  i.e.  the  presence  of  a  maximum  in  the  time  to 
explosion,  which  delineates  two  regimes.  The  peak  in  texp  for  MP5  in  Figure  6.7 
occurs  at  approximately  Ax  =  7  fim,  whereas  the  critical  Ax  value  is  slightly 

lower,  at  approximately  5.5  /im,  for  AW5,  and  event  lowere  for  MW5  at  Ax  ~  Afim. 
For  grid  resolutions  Ax  below  the  critical  value,  the  numerical  simulation  is  in  a 
“convectively”  dominant  regime,  where  the  combination  of  the  numerical  scheme 
and  fine  grid  resolution  is  sufficient  to  effectively  mitigate  the  effects  of  numerical 
diffusion  in  the  detonation  formation.  AW5  &  MW5  produce  similar  values  of  texp 
for  Ax  <  3  /im,  but  for  this  complex  kinetics  scheme,  time  convergence  for  texp 
may  not  be  reached  except  for  Ax  <  0.5  fim,  consistent  with  findings  by  Powers 
&  Paolucci[72],  In  the  case  of  AW5,  grid  convergence  is  observed  at  Ax  ~  0.5 fim. 

125 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


For  grid  resolutions  greater  than  the  critical  value  of  Ax  in  Figure  6.7,  we  enter  the 
numerically  dissipative  regime,  where  coupling  of  the  fluid  mechanics  and  kinetics 
is  enhanced  due  to  numerical  diffusion  of  temperature  and  chemical  concentrations. 
The  value  of  texp  thus  decreases  with  increasing  Ax  in  this  numerically  dissipative 
regime. 

To  obtain  results  that  are  truly  insensitive  to  the  numerical  effects,  one  requires 
Ax  — *  0,  since  the  two  methods  converge  in  that  limit  to  a  single  value  for  texp,  but 
of  course  this  is  a  practical  impossibility.  The  results  in  Figure  6.7  suggest  that, 
to  be  able  to  reasonably  resolve  detonation  propagation  and  shock-flame  coupling 
dynamics,  a  grid  spacing  of  Ax  =  2.5  /irn  or  smaller  may  produce  acceptable  ac¬ 
curacy,  given  less  than  a  5%  difference  between  the  two  different  schemes.  Yet  the 
time  to  re-explosion  is  but  one  parameter  that  is  affected  by  grid  resolution,  and 
calculations  of  the  shock-flame  instabilities  at  smaller  grid  spacings  are  needed  to  be 
able  to  explore  other  quantitative  features. 

A  simple  fast  Fourier  transform  (FFT)  can  be  used  to  find  the  spectral  content 
of  the  two  instability  modes,  HF  and  HA,  which  are  observed  in  Figures  6.2  -  6.4. 
Figure  6.8  shows  FFT  results  for  the  MP5  scheme  and  grid  sizes  Ax  =  2.5,  1.5,  and 
1.0  /irn.  There  is  remarkable  consistency  in  the  dominant  frequencies  here;  for  both 
HA  instabilities  near  a  frequency  of  approximately  0.35  MHz  and  HF  instabilities 
near  2.3  MHz,  there  is  relatively  little  difference  in  results  for  Ax  <  2.5/im.  This 
observation  suggests  that  the  spectral  content  of  the  resulting  instability  is  relatively 
insensitive  to  the  grid  spacing,  as  long  as  Ax  lies  below  the  critical  value  for  the 
start  of  numerical  diffusion.  At  high  frequencies,  above  4  MHz,  there  is  some  grid 
dependency;  in  fact,  frequencies  above  4  MHz  are  not  seen  for  Ax  =  1.5  fim.  The 

spectral  content  around  4.5  MHz  is  likely  to  be  a  harmonic  of  the  strong  2.3  MHz 
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signal.  Note  that  since  the  instabilities  develop  for  a  finite  time  only,  the  sampling 
statistics  of  the  FFT  are  limited.  The  use  of  a  wavelet  decomposition  did  not  provide 
improvements  in  the  signal-to-noise  ratio. 

These  results  confirm  earlier  complex  kinetics  findings  [16]  which  indicate  that 
a  detonation  near  the  CJ  limit  has  two  physically  distinct  instability  modes.  The 
high  frequency  mode  always  appears  first  and  marks  the  transition  from  a  ‘stable’ 
CJ  detonation,  where  a  low  frequency  mode  appears  later  in  time.  The  overdriven 
detonation  simulations  by  Leung  et  al.  [IT]  with  a  two-step  reaction  mechanism  also 
demonstrate  multiple  instability  modes,  but  with  distinct  differences  from  those  seen 
here,  as  noted  above. 

Our  observed  fluctuations  in  key  properties  (e.g.,  in  species  concentration,  tem¬ 
perature,  and  pressure)  of  the  fluid  within  the  induction  zone  are  described  by  Oran 
&  Boris[73]  as  ‘hot  spots’.  The  present  study  with  complex  reaction  kinetics  shows 
that  these  ‘hot  spots’  contribute  to  an  initial  stage  of  the  flame  dynamics.  In  this 
regime,  the  induction  length  is  very  small  (£  «  £cj),  and  acoustic  waves  gener¬ 
ated  by  the  perturbed  chemistry  are  rapidly  transmitted  to  the  shock,  i.e.  leading 
to  high-frequency  modes.  Because  there  is  a  very  limited  amount  of  fluid  that  can 
participate  in  the  fluctuation  of  the  heat  release,  only  low-amplitude  perturbations 
of  the  CJ  peak  pressure  appear.  As  these  acoustic  waves  reach  the  leading  shock  and 
strengthen  it,  their  frequency  can  be  measured  as  that  of  the  fluctuations  of  the  peak 
pressure.  Eventually  the  average  induction  length  continues  to  increase  and  the  sec¬ 
ond  mode  appears,  which  directly  couples  the  flame  speed  with  the  shock,  resulting 
in  fluctuations  with  lower  frequency  but  much  higher  amplitude.  This  interpretation 
of  our  observations  will  be  elaborated  upon  and  verified  in  the  next  section. 
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6.2  Simplified  Model 


To  better  understand  and  interpret  the  coupling  between  reactive  and  fluid  me- 
chanical/acoustic  phenomena  that  generates  our  observed  results,  a  model  for  the 
induction  zone  may  be  constructed  and  explored.  This  model  is  composed  of  a  lead¬ 
ing  shock,  a  heated,  post-shock  medium  (fluid),  and  a  flame  front,  all  of  which  are 
illustrated  in  Figure  6.9.  A  single  period  of  the  detonation  oscillation  can  be  de¬ 
scribed,  in  the  reference  frame  of  the  shock,  in  a  manner  similar  to  that  of  McVey  & 
Toong[74],  as  follows.  Fluctuations  at  the  flame  front  create  an  acoustic  (pressure) 
disturbance,  which  travels  at  the  acoustic  wave  speed,  \ac,  through  the  induction 
zone  until  it  reaches  the  leading  shock;  this  process  occurs  between  reference  times 
ta  and  tb .  Upon  contact,  the  pressure  fluctuation  carried  by  the  acoustic  wave  will 
accelerate  the  shock  and  alter  the  post-shock  conditions,  thus  creating  an  entropy 
disturbance  (temperature  fluctuation).  This  entropy  disturbance  will  propagate  back 
into  the  induction  zone  at  the  entropy  wave  speed,  Aen,  toward  the  flame  front,  this 
process  occuring  between  times  tb  and  tc.  A  resonant  condition  is  achieved  when, 
upon  contact  with  the  flame,  the  entropy  wave  creates  a  new  acoustic  disturbance 
in  the  flame,  and  the  cycle  repeats.  Figure  6.9  illustrates  this  phenomenon,  with 
relations  for  the  entropy  and  acoustic  wave  speeds  as  follows,  respectively: 


.  .  dx 

Kn(x,i)  =  ~Jt 


.  ,  .  dx 

A„(x,«)=  - 


=  u2(x,t) 

en 

c(x,t)  -  u2{x,t) 


(6.1) 

(6.2) 


Here  u(x,t )  is  the  fluid  velocity,  c(x,t)  is  local  speed  of  sound,  D{t)  is  the  detona¬ 
tion  velocity,  and  u2(x,t )  =  | u(x,t)  —  D(t) |  is  the  post-shock  fluid  velocity  in  the 
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detonation  reference  frame.  From  these  wave  speeds  the  period  of  the  cycle,  r,  may 
be  expressed  by 


where  xs 
flame. 


r  = 


lxf 


-^ac("D  £) 


dx 


+ 


era  (•£ ,  t) 


dx 


(6,3) 


(t  —  to)  ■  D(t )  is  the  position  of  the  shock  and  Xf  is  the  position  of  the 


At  a  zeroth-order  approximation,  the  fluid  properties  in  the  induction  region, 
Q2(x,t),  are  assumed  to  weakly  vary  with  time  for  a  given  half  cycle,  ~  0  & 

—  0-  This  also  implies  that  the  relative  positions  of  the  flame  and  shock  front 
are  approximately  constant.  From  this  approximation,  the  period  can  be  determined 


C-ab  T  ^ ab  Dab  ^ be  D^c 


(6.4) 


where  i  is  the  period-averaged  induction  length,  ab  is  the  fluid  state  at  the  acoustic 
wave  half-cycle,  be  is  the  fluid  state  during  the  entropy  wave  half  cycle,  and  ua,  ca, 
and  Da  are  the  the  fluid  speed  in  the  detonation  reference  frame,  speed  of  sound, 
and  average  detonation  speed,  respectively,  for  half-cycle  a.  The  model  for  acoustic 
and  entropy  half  cycles  are  illustrated  in  Figure  6.10,  corresponding  to  observed 
oscillations  as  indicated  in  the  inset.  From  the  period  of  the  combined  cycles,  the 
frequency  in  oscillations  of  the  peak  pressure  trace  is  /  =  r^1. 


Data  may  be  extracted  from  the  full  kinetics  results  at  different  peak  pressure 
cycles  from  the  high  frequency  as  well  as  high  amplitude  regime,  and  compared 
with  the  simplified  model  expressed  in  Equation  6.4.  Figure  6.12  illustrates  the 
evolution  of  the  induction  zone  temperature  profile  in  the  detonation  reference  frame 

for  a  given  period  of  the  high  amplitude  mode,  with  data  for  the  acoustic  wave 
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propagation  in  (a)  and  for  the  entropy  wave  in  (b).  The  density  distributions  for 
the  corresponding  HA  mode  are  shown  in  Figure  6.14.  Using  the  induction  zone 
data  and  the  period  from  Eq.  (6.4),  the  frequency  /  m  310  kHz  may  estimated  for 
the  HA  mode,  which  is  in  excellent  agreement  with  that  obtained  from  the  spectral 
analysis  using  the  full  simulation,  310  ±  40  kHz,  shown  in  Figure  6.8.  Performing 
the  same  analysis  on  the  high  frequency  (HF)  mode,  illustrated  in  the  temperature 
profiles  in  Figure  6.11  and  the  density  profiles  in  Figure  6.13,  the  frequency  /  ~ 
2.08  MHz  may  estimated,  which  is  also  in  good  agreement  with  that  extracted 
from  the  spectral  analysis  (2.29  ±  0.4  Mhz )  in  Figure  6.8.  Hence  this  simple  model 
appears  to  capture  reasonably  well  the  global  processes  that  lead  to  the  high-  and 
low-frequency  detonation  instability  modes. 

We  should  point  out  that  this  simple  model  has  been  succesfully  applied  to  an¬ 
other  case  of  shock-induced  instability.  For  strong,  ionizing  shocks  in  a  noble  gas 
(specifically,  argon),  a  similar  structure  of  shock,  induction  zone  and  reaction  front 
can  be  observed  [75,  76,  77].  The  instability  in  this  case  has  lower  amplitudes,  due  to 
the  absence  of  exo-thermic  reactions,  but  is  also  well  explained  by  a  resonant  coupling 
between  the  shock  front  and  reaction  zone  through  the  transmission  and  reflection  of 
entropy  and  acoustic  waves  [75].  In  fact,  the  2D  structure  of  that  flow  could  be  found 
to  have  features  remarkably  similar  to  detonations  including,  for  example,  artihcal 
soot  patterns  [77].  This  indicates  that  the  model,  despite  its  simplicity,  provides  a 
good  insight  into  the  true  physical  mechanisms  involved  and  may  have  relevance  to 
a  universal  class  of  instabilities  in  reactive  shock  systems. 
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6.3  Discussion 


The  ‘hot  spot’  which  appears  in  the  high  frequency  mode  temperature  profiles  (Figure 
6.11)  is  of  particular  interest  in  interpreting  the  differences  between  modes.  This  hot 
spot  burns  only  a  fraction  of  the  overall  mixture  and  the  heat  release  is  not  sufficient 
to  significantly  alter  the  characteristic  speeds  of  the  flow  (especially  since  the  speed 
of  sound  only  varies  with  the  square-root  of  the  temperature).  Examination  of  the 
temperature  profiles  in  Figure  6.11,  in  particular  panel  (a),  shows  that  the  reaction 
zone  fluctuates  from  a  ‘hot-spot7  position  without  significant  change  in  location;  if 
at  all,  the  perturbed  region  of  accelerated  burning  is  actually  convected  downstream. 
This  pre-ignition  effect  does  not  allow  the  flame  to  accelerate,  and  in  the  case  of 
the  high  frequency  mode,  the  detonation  is  still  slightly  overdriven.  The  resonance 
between  the  perturbation  of  the  chemical  rates  at  the  hot  spot  and  the  shock  front 
still  remains,  though,  and  is  the  basis  for  the  observed  oscillation  pattern,  as  indicated 
by  the  agreement  between  computed  and  measured  frequencies. 

By  contrast  in  the  high  amplitude  mode’s  temperature  profiles,  seen  in  Figure 
6.12,  there  is  no  observed  ‘hot  spot7.  More  likely,  the  perturbation  has  moved  to  the 
flame  region  and  any  acceleration  of  the  chemical  rates  in  that  region  can  enhance 
the  rate  of  heat  release  much  more  significantly.  This  allows  what  is  presumably  a 
SWACER-like  [71]  mechanism  to  govern  this  high  amplitude  regime,  starting  from 
a  flame  at  the  furthest  distance  from  the  shock,  £ma3.,  and  accelerating  towards  the 
shock  as  it  burns  the  fluid  in  the  induction  zone,  releasing  large  amounts  of  energy 
up  to  the  point  of  contact  with  the  shock.  This  can  be  seen  in  Figure  6.12(a),  where 
the  temperature  profiles  for  the  HA  mode  clearly  exhibit  a  dramatic  reduction  of  the 
distance  between  the  shock  and  flame  fronts  in  time.  Thus,  the  HA  mode  is  very 
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much  a  SWACER-like  mechanism,  while  the  HF  mode  is  not.  This  large  change  in 
shock-flame  distance  is  not  seen  in  the  two-step  model  results  of  Leung  et  al.  [17], 
thus  highlighting  another  important  feature  of  the  true  kinetics. 

From  the  density  profiles  illustrated  in  Figures  6.13  and  6.14,  it  becomes  quite 
apparent,  even  in  the  high-frequency  regime,  that  there  is  a  large  variation  of  density 
within  the  induction  zone  throughout  a  given  cycle.  Neither  single-step  Arrhenius 
kinetics  [12,  13,  14,  15]  nor  a  two-step  reaction  model  [17]  account  for  these  large 
density  variations;  their  ZND  approximation  is  only  valid  for  a  brief  portion  of  the 
cycle  in  the  HA  regime.  Therefore,  more  complex  reaction  kinetics  than  have  been 
incorporated  in  the  past  should  be  utilized  to  capture  the  full  quantitative  features  of 
the  evolving  detonation  instabilities.  The  simple  two-wave  resonant  model  explored 
in  Section  6.2  also  depends  on  approximately  constant  flow  properties  within  the  in¬ 
duction  zone,  and  thus  the  same  limitations  apply  to  this  model,  except  that  we  have 
separated  the  cycle  into  two  sections,  each  with  different  average  flow  properties  and 
induction  lengths.  Hence  a  reasonable  representation  of  global  dynamical  character 
is  achieved  with  the  simple  model. 

As  the  detonation  relaxes  toward  the  CJ  condition  prior  to  the  onset  of  the 
instabilities,  the  post  shock  conditions  can  be  used  to  determine  the  chemical  time 
scale  corresponding  to  the  induction  zone.  For  the  present  simulations  this  time  scale 
is  approximately  Thydr  =  300  ±  10ns.  Using  the  same  post  shock  conditions  in  a  zero¬ 
dimensional  reactive  simulation  produces  time  scale  Tchem  =  215  ±  5 ns.  This  level  of 
disparity  is  to  be  expected,  because  the  zero- dimensional  calculations  are  performed 
at  constant  volume,  and  as  the  system  approaches  the  peak  of  the  H  concentration 
(the  criterion  used  to  define  the  flame  center),  the  energy  liberated  remains  confined 

and  leads  to  a  more  rapid  rate  of  reaction.  Using  the  Thydr  value  as  a  normalization 
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factor,  the  two  modes  observed  in  the  present  simulations,  HF  and  HA,  correspond 
to  normalized  periods  of  1.6  and  10.75,  respectively.  This  is  within  the  range  of 
instabilities  observed  by  Leung  et  al.  [17] ;  a  more  exact  correspondence  is  difficult  to 
obtain,  since  several  parameters  are  being  varied  in  their  two-step  model,  with  no 
direct  relation  to  the  actual  chemical  system.  As  noted  in  Section  6.1,  in  the  present 
studies  we  do  not  observe  the  so-called  “very-high  frequency”  mode  seen  by  Leung  et 
al.[17],  except  for  very  low  amplitude  oscillations  that  are  shown  to  be  a  numerical 
artifact.  In  fact,  a  physical  instability  at  very  high  frequency  would  correspond  to  an 
oscillation  period  that  is  less  than  the  induction  delay,  implying  that  the  “hot  spot” 
be  located  very  close  to  the  shock.  In  that  case,  the  extent  of  this  perturbation  (i.e., 
the  reaction  time)  would  need  to  be  much  smaller  than  the  induction  period.  With 
realistic  chemistry,  this  may  not  be  possible.  Thus,  it  is  conjectured  here  that  the 
very  high  frequency  mode  may  be  an  artifice  of  the  two-step  model. 
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(a)  PSpark  =  20  atm  with  0.25  cm  spark  length  where  detonation  is  not  achieved 


x  [cm] 

(b)  PSpark  =  40  atm  with  0.25  cm  spark  length  where  detonation  is  achieved 


Figure  6.1:  Pressure  contours  on  an  x-t  diagram  for  a  spark  ignited  H2-A1T  mixture 
with  Ax  =  50 fim,  computed  using  the  MP5  scheme. 
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Figure  6.2:  Peak  pressure-time  history  of  a  spark-ignited  i^-air  mixture  simulated 
with  two  different  grid  cell  sizes  Ax  using  the  MP5  scheme.  The  time  to  re-explosion, 
texp ,  high  amplitude  mode,  HA,  and  high  frequency  mode,  HF,  are  illustrated. 
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(a)  Ax  =  12.5  [im 


(b)  Ax  =  2.5  fj,m 

Figure  6.3:  High  frequency  portion  of  peak  pressure  time  history  as  in  Fig.  6.2, 
simulated  with  two  different  grid  cell  sizes  Ax. 
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(a)  Ax  =  12.5  /jm 


(b)  Ax  =  2.5  fj,m 

Figure  6.4:  High  amplitude  portion  of  peak  pressure  time  history  as  in  Fig.  6.2, 
simulated  with  two  different  grid  cell  sizes  Ax. 
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Figure  6.5:  Early  portion  of  peak  pressure  time  history  after  re-explosion,  simulated 
with  two  different  grid  cell  sizes  Ax. 
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(a)  Induction  length  vs.  time 


(b)  Induction  lengh  vs.  time  (zoom) 


Figure  6.6:  Early  portion  of  induction  length  time  history  after  re-explosion,  simu¬ 
lated  with  two  different  grid  cell  sizes  Ax. 


139 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


Figure  6.7:  Time  to  re-explosion  as  a  function  of  grid  resolution  size  Ax  for  the  MP5, 
AW5,  and  MW5  schemes. 
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Scaled 


Figure  6.8:  Fast  Fourier  transform  of  peak  pressure  traces  for  the  unstable  detonation 
with  3  different  grid  sizes  Ax,  simulated  using  the  MP5  scheme. 
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Figure  6.9:  A  simple  model  involving  the  shock,  flame,  induction  zone,  and  trans¬ 
mission  of  entropy  and  acoustic  waves  to  represent  shock-flame  coupling. 


Figure  6.10:  Simplified  model  of  the  peak  pressure  cycle  (left)  used  to  represent  the 
numerically  observed  pulsations  in  peak  pressure  (right  inset). 
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3000 


x  [mm] 

(a)  acoustic  wave  cycle 


(b)  entropy  wave  cycle 

Figure  6.11:  Temperature  distribution  in  shock  reference  frame  at  different  times 
within:  (a)  the  acoustic  wave  cycle  (times  a  =  29.1  /is  to  b  =  29.3  fis)  and  (b)  the 
entropy  wave  cycle  (times  b  =  29.3  /is  to  c  =  29.5  /is).  Data  are  extracted  from  high 
frequency  (HF)  mode  results. 
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x  [mm] 

(a)  acoustic  wave  cycle 


(b)  entropy  wave  cycle 


Figure  6.12:  Temperature  distribution  in  shock  reference  frame  at  different  times 
within:  (a)  the  acoustic  wave  cycle  (times  a  =  57.1  /is  to  b  =  58.3  {is)  and  (b)  the 
entropy  wave  cycle  (times  b  =  58.3  {is  to  c  =  60.7  {is )  Data  are  extracted  from  high 
amplitude  (HA)  mode  results. 

144 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


5.5 


x  [mm] 

(a)  acoustic  wave  cycle 


(b)  entropy  wave  cycle 

Figure  6.13:  Density  distribution  in  shock  reference  frame  at  different  times  within: 
(a)  the  acoustic  wave  cycle  (times  a  =  29.1  (is  to  b  =  29.3  (is)  and  (b)  the  entropy 
wave  cycle  (times  b  =  29.3  (is  to  c  =  29.5  (is).  Data  are  extracted  from  high  frequency 
(HF)  mode  results. 
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6 


(a)  acoustic  wave  cycle 


(b)  entropy  wave  cycle 

Figure  6.14:  Density  distribution  in  shock  reference  frame  at  different  times  within: 
(a)  the  acoustic  wave  cycle  (times  a  =  57.1  /is  to  b  =  58.3  /is)  and  (b)  the  entropy 
wave  cycle  (times  b  =  58.3  /is  to  c  =  60.7  /is).  Data  are  extracted  from  high 
amplitude  (HA)  mode  results. 
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CHAPTER  7 


Magnetic  Field  and  Detonation  Interactions 


In  evaluating  the  feasibility  of  the  PDRIME  configurations  we  must  examine  how  a 
magnetic  held  will  interact  with  a  conducting,  reacting  how  such  as  the  detonation 
as  well  as  the  amount  of  conductivity  required  in  the  fluid  for  sufficient  interaction 
of  the  magnetic  held  and  detonation.  In  Chapter  5,  various  PDRIME  configurations 
have  been  evaluated  using  simplified  combustion  kinetics  and  the  exploration  of  the 
effect  of  MHD  to  accelerate  the  reactive  how  in  the  bypass  section  or  magnetic 
chamber  piston  concept  (CP)  or  decelerate  the  how  in  the  nozzle  via  nozzle  generator 
concept(NG).  In  these  prior  calculations,  the  ability  of  MHD  to  affect  combustion 
was  modeled  globally,  but  not  in  detail.  The  goal  of  the  studies  in  this  chapter  is 
to  determine  if  such  MHD-based  acceleration/deceleration  is  possible  for  detonation 
phenomena. 

7.1  Detonation  Instabilities  with  Applied  Magnetic  Fields 

In  Section  6.1,  the  stability  of  an  unsupported  detonation  was  evaluated.  For  the 
next  simulations,  we  began  with  the  same  initial  conditions  used  in  the  previous 
chapter  and  seeded  the  fluid  with  different  amounts  of  cesium  (1%,  5%,  and  10%  by 
mole).  Because  of  the  large  molecular  weight  of  cesium  (f^  ~  70,  ~  5),  adding 
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upwards  of  10%  Cs  can  significantly  affect  the  dynamics  of  the  detonation  (e.g., 
CJ  detonation  speed  Dqj ,  flame  temperature  Tf,  post-leading  shock  pressure  and 
temperature  PvN ,TvN ,  etc.)  even  without  the  presence  of  an  applied  magnetic  held. 
Figure  7.1  illustrates  the  peak  pressure  traces  of  unmagnetized  spark-ignited  H2— air 
detonations  with  the  addition  of  0,  1,  5,  &  10  %  Cs.  A  grid  resolution  of  5/im  and  the 
MP5  scheme  were  used  to  perform  the  ID  detonation  simulations  in  the  remainder 
of  this  section.  Figure  7.1(a),  in  the  absence  of  Cs  addition,  shows  the  initiation 
of  HF  and  HA  instability  modes,  as  observed  in  Figures  6.2(a)  &  (b)  for  different 
grid  cell  sizes.  With  a  grid  resolution  of  5  /im  it  is  possible  to  capture  the  essential 
dynamics  of  the  detonation,  within  the  “convection”  dominated  regime  (see  Figure 
6.7),  although  parameters  such  as  time  to  re-explosion  may  not  be  as  accurately 
determined.  Nevertheless,  5  /im  resolution  is  sufficient  to  enable  computation  of 
additional  ionization  processes  of  Cs  and  their  effect  on  instabilities. 

From  Figure  7.1  (a)-(d),  one  will  notice  the  oscillations  become  less  erratic  and,  at 
5%  and  10%  molar  addition,  approaching  consistent  HA  behavior,  as  more  cesium  is 
added.  Extinction,  where  the  peak  pressure  decays,  is  significantly  delayed  for  higher 
Cs  concentrations.  This  trend  is  consistent  with  the  observations  of  Radulescu  et 
al.  [19]  where  argon  was  incrementally  added  to  acetylene-oxygen  mixture  to  stabilize 
the  detonation.  Heavy  argon  dilution  in  the  mixture  led  to  large-frequency,  small 
amplitude  regular  oscillations  of  the  shock  front  pressure.  Radulescu  et  al.  attributed 
the  stabilizing  effect  of  the  diluent  to  the  lower  temperature  in  the  reaction  zone  which 
leads  to  slower  exothermic  reaction  rates.  For  the  present  calculations,  regularization 
of  both  the  HF  and  HA  modes  with  Cs  addition  could  result  from  the  lowering  of 
the  reaction  front  temperature.  The  addition  of  Cs  is  also  observed  to  reduce  the 
effective  ID  detonation  speed,  as  also  determined  by  the  theoretical  reduction  in  the 
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CJ  speed  (see  Table  7.1).  The  results  shown  in  Figure  7.1  will  serve  as  the  “baseline” 
case  for  the  ID  dynamics  as  affected  by  an  applied  magnetic  field. 

Now  that  the  “baseline”  behavior  has  been  established  for  cesium-seeded  det¬ 
onations,  let  us  now  examine  how  the  behavior  of  these  instabilities  change  when 
one  gradually  increases  the  strength  of  an  applied  magnetic  field.  A  fixed  magnetic 
field  is  used,  per  Equation  2.17,  with  a  loading  factor,  K  =  0  (no  applied  electric 
field).  The  magnetic  field  configuration  is  illustrated  in  Figure  7.2  where  a  transverse 
magnetic  field  will  be  applied  with  strengths  ranging  from  0  to  8  Tesla.  A  spark- 
ignited  detonation  will  propagate  into  the  region  of  this  magnetic  field,  as  shown. 
Figures  7.3,  7.4,  and  7.5  show  the  peak  pressure  traces  of  detonations  seeded  with 
1%,  5%,  and  10%  Cs,  respectively,  under  the  influence  of  various  applied  magnetic 
field  strengths.  Without  an  applied  magnetic  field,  the  peak  pressure  trace  of  a  det¬ 
onation  with  a  1%  Cs  has  an  irregular  periodicity  (Figures  7.1(b)  &  7.3(a)).  As  the 
magnetic  field  is  strengthened  (see  Figures  7.3(b)-(f))  the  detonation  peak  pressure 
trace  shows  more  prolonged  erratic  behavior,  with  no  clear  trend  in  the  time  at  which 
extinction  (pressure  decay)  takes  place.  But  Figure  7.4  shows  that  as  the  magnetic 
field  strength  is  increased,  the  oscillations  of  the  peak  pressure  trace  are  driven  to 
become  less  regular  at  earlier  times  and  extinguish  sooner.  Similarly  findings  are 
observed  in  Figure  7.5.  There  are  a  few  qualitative  assessments  that  can  be  made 
from  these  results.  Figures  7.3(a),  7.4(d),  &  7.5(e)  illustrate  peak  pressure  traces 
of  detonations  seeded  with  1%,  5%,  and  10%  cesium  under  applied  field  strengths 
of  Bz  =  0,  6,  and  7  Tesla,  respectively.  These  peak  pressure  traces  demonstrate  a 
general  trend.  It  would  appear  as  though  the  stabilizing  effect  of  the  presence  of 
diluent,  cesium,  becomes  less  effective  when  the  fluid  at  the  flame  has  a  greater  con¬ 
ductivity  and  thus  introduces  an  additional  scale  to  the  flame-shock  dynamics.  One 
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could  attempt  to  find  the  modes  that  exist  within  the  peak  pressure  trace  to  give 
some  quantitative  insight  into  the  effects  of  these  fields,  but  the  data  would  prove 
unreliable  due  to  the  small  time  frame  during  which  these  large  oscillations  exist. 
The  fact  that  the  application  of  the  magnetic  held  can  degrade  the  cyclic  stability 
of  the  ID  detonation  at  earlier  times  does  indicate  the  ability  of  the  B  held  to  alter 
the  hame-detonation  coupling  as  well  as  the  combustion  process  itself. 

7.2  Detonation  Instabilities  with  MHD 

The  behavior  of  a  detonation  with  an  applied  held  using  the  MHD  accelerator  con¬ 
figuration  will  now  be  investigated.  The  study  will  begin  with  a  hxed  magnetic  held, 
Equation  2.17,  and  the  loading  factors  recommended  by  Cambier[4],  Ky  =  1.5  and 
0.5  (i.e.,  the  electric  held  scaled  with  ux  x  Bz ),  for  the  accelerator  and  generator, 
respectively.  The  remainder  of  the  current  detonation  configuration  is  as  specihed  in 
Section  7.1. 

The  previous  section  demonstrated  that  the  H2-& ir,  ID  spark-induced  detonation 
under  investigation  is  unstable  under  normal  operating  conditions  with  or  without  an 
applied  magnetic  held.  Figures  7.6-7. 8  show  the  peak  pressure  traces  of  detonation 
seeded  with  1%,  5%,  and  10%  Cs,  respectively,  with  an  electric  held  in  the  acceler¬ 
ator  configuration  (Ky  =  1.5)  under  the  influence  of  various  applied  magnetic  held 
strengths.  In  the  1%  Cs  test  case,  it  becomes  quite  clear,  when  comparing  Figures 
7.6(b)-(f)  with  the  unmagnetized  case  of  Figure  7.6(a)  that  the  MHD  acceleration 
has  the  ability  to  regularize  and  sustain  the  oscillating  detonation;  the  detonation 
does  not  extinguish,  for  the  time  period  shown,  as  it  does  in  Figure  7.6(a).  When 
the  amount  of  cesium  is  increased  to  5%,  as  the  magnetic  held  is  increased  as  shown 
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in  Figures  7.7(a)-(f)  the  detonation  is  cyclically  stabilized  and  the  peak  pressure 
amplitude  greatly  reduced,  by  around  30%.  When  the  amount  of  cesium  is  increased 
to  10%,  as  the  magnetic  field  is  increased  from  Bz  =  0 T,  Figure  7.8(a),  to  Bz  =  8 T, 
Figure  7.8(f),  the  detonation  goes  from  a  cyclically  stable  galloping  detonation  to  a 
nearly  stable  CJ  detonation,  with  a  slightly  delayed  onset  of  the  initial  instabilities. 

Because  of  the  strong  effect  the  magnetic  field  had  on  the  detonation  stability 
and  sustainment,  in  both  the  accelerator  and  generator  configurations,  we  now  ex¬ 
amine  these  features  side-by-side.  Figures  7.9(a) (c)(e)  illustrate  the  time  dependent 
peak  pressure,  induction  length,  and  detonation  velocity,  respectively,  for  a  10  %  Cs 
in  the  accelerator  configuration  with  applied  magnetic  field  strengths  of  0,  3,  and  8 
T.  Figures  7.9(b) (d)(f)  show  the  same  thing  for  the  generator  configuration.  Fig¬ 
ure  7.9(a)  &  (c)  show  a  significant  decrease  in  amplitude  of  the  peak  pressure  and 
induction  length,  respectively,  as  the  magnetic  field  strength  is  increased  for  an  accel¬ 
erator,  while  figure  7.9(e)  shows  the  detonation  speed  stabilizes  near  but  above  the 
Chapman- Jouguet  velocity,  with  little  propensity  for  extinction.  In  contrast,  Figures 
7.9(b)  &  (f)  for  the  generator  configuration  show  that  peak  pressure  and  detonation 
speed  become  erratic,  with  increasing  magnetic  field  strength,  and  Figure  7.9(d) 
shows  that  the  induction  length  oscillations  becomes  unbounded,  leading  to  extinc¬ 
tion  for  the  8  Tesla  case.  This  contrast  is  quite  striking,  especially  in  the  generator’s 
influence  on  extinguishing  the  detonation.  There  appears  to  be  some  acceleration  of 
the  detonation,  to  a  speed  greater  than  CJ  for  the  accelerator  configuration  (Figure 
7.9(e),  8T),  but  overall  there  is  no  significant  difference  in  the  detonation  velocities 
between  accelerator  and  generator  configurations. 
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7.3  2D  Cellular  Detonation  in  an  Applied  Magnetic  Field 


While  the  effects  of  MHD  on  the  stability  of  the  ID  detonation  are  interesting,  the 
more  realistic  flow  is  that  of  a  two-dimensional  detonation  structure.  Before  the 
effects  of  an  applied  magnetic  field  on  a  two-dimensional  detonation  are  examined, 
however,  one  must  first  examine  the  effects  of  seeding  cesium  on  the  evolving  cellular 
detonation.  The  detonation  front  as  well  as  the  underlying  detonation  structure  may 
be  found  in  the  contours  of  the  maximum  pressure,  Pmax ,  at  each  point  in  space. 
These  maximum  pressure  traces  have  been  observed  in  experiments  over  decades  in 
smoke-foil  records[78].  The  unmagnetized  2D  detonation  with  various  amounts  of 
cesium  will  be  examined  first. 

The  2D  detonation  initiated  with  a  computational  spark  is  illustrated  in  Figure 
7.10.  A  grid  resolution  of  Ax  =  Ay  =  50 ym  and  the  MW5  scheme  are  used  to 
perform  the  2D  detonation  simulations  in  the  remainder  of  this  chapter;  much  finer 
grid  resolution  becomes  prohibitively  expensive,  yet  this  is  found  to  be  sufficient  for 
the  study  of  the  overall  detonation  structure.  Figure  7.11  shows  that  evolution  of 
an  H2-a.i1  detonation  front  (Schlieren-like  plot  and  smoke-foil  record)  without  the 
presence  of  an  applied  magnetic  field  and  without  Cs  injection.  From  the  series  of 
figures,  the  Mach  stem,  transverse  shocks,  and  incidents  shocks  are  well  resolved. 
Figure  7.12  shows  the  detonation  front  and  smoke-foil  record  after  75 ys  of  the  same 
detonation.  At  this  resolution,  it  can  be  clearly  seen  in  the  smoke-foil  record,  Figure 
7.12(b),  that  cellular  detonation  is  achieved,  and  is  consistent  with  the  established 
literature  for  cellular  detonations [78,  79,  80,  81].  The  simulated  detonation  cellular 
length  and  width  for  this  configuration  are  A l  =  2.7±0.1  mm  and  Aw  =  1.67±0.1mm, 
respectively.  It  can  also  be  demonstrated  that  an  unmagnetized  H2-air-Cs  mixture 
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seeded  with  cesium  ranging  from  1-10%  can  achieve  the  same  cellular  detonation  and 
detonation  front  resolution,  shown  in  Figures  7.13  -7.15.  The  decreased  detonation 
velocity  as  cesium  is  increased  from  1%  to  10%  is  consistent  with  the  ID  simulations 
as  well  as  the  theoretical  CJ  velocity,  which  are  tabulated  in  Table  7.1. 

In  the  previous  sections,  7.1  and  7.2,  we  demonstrated  that  an  applied  field  can 
alter  the  cyclic  stability  of  a  ID  detonation.  Now  we  will  investigate  if  an  applied 
field  can  affect  the  regularity  or  velocity  of  a  2D  cellular  detonation.  Due  to  the  large 
spatial  variation  of  the  longitudinal  velocity,  ux ,  the  static  electric  field  should  no 
longer  be  optimally  tailored  to  ux  and  Bz.  In  order  to  more  accurately  account  for 
the  static  electric  fields  contribution  to  the  MHD  forces,  while  still  keeping  the  load¬ 
ing  factor  concept  for  the  accelerator  and  generator  configurations,  the  static  electric 
field  is  prescribed  as  Ey  =  KyU Bz,  where  U  is  the  mean  velocity,  and  will  not  evolve 
with  the  flow.  The  optimal  condition  U  ~  1000  m/s  is  assumed.  As  performed 
with  the  one- dimensional  test  cases,  2D  detonations  in  the  both  the  accelerator  and 
generator  configurations  were  computed  with  varying  amounts  of  cesium(l,  5,  and 
10%)  using  the  highest  magnetic  field  strength  used  previously  in  the  ID  detona¬ 
tion  simulations,  Bz  =  8 T.  From  the  X-t  plots  of  the  centerline  of  the  2D  shock 
front,  illustrated  in  Figures  7.16-7.18,  there  is  no  noticeable/significant  altering  of 
the  detonation  velocity  in  either  generator  or  accelerator  modes,  as  compared  with 
the  detonation  in  the  absence  of  MHD.  The  peak  pressure  evolution  in  Figures  7.19- 
7.21  similarly  do  not  show  any  noticeable  change  in  the  cyclic  stability  of  these 
detonations.  The  differences  in  the  two  configurations  can  be  seen  in  the  centerline 
pressure  and  conductivity  profiles,  where  the  10  %  Cs  and  B,  =  8 T  case  at  t  =  75{is 
case  is  illustrated  in  Figure  7.22.  These  differences  in  the  profiles  are  seen  further 
downstream,  however,  and  do  not  play  a  role  in  the  detonation  dynamics.  Hence  even 
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at  a  very  high  magnetic  field  strength,  8T,  and  for  10%  Cs,  there  does  not  appear 
to  be  a  significant  influence  of  the  MHD  on  detonation  dynamics.  In  an  attempt  to 
see  if  the  two-dimensional  detonation  can  be  altered  by  MHD  in  some  meaningful 
way,  but  under  different  conditions,  two  alternative  approaches  were  taken.  The  first 
was  to  artificially  enhance  the  kinetics  of  the  cesium,  (the  original  formulation  of 
which  is  described  in  Appendix  A).  The  second  approach  drastically  increased  the 
strengths  of  the  applied  magnetic  and  electric  fields.  For  both  cases,  10  %  Cs  was 
seeded  because  it  allowed  for  more  optimal  levels  of  ionization  as  compared  to  1  and 
5%  Cs  addition. 

7.3.1  Enhanced  Kinetics 

In  the  enhanced  kinetics  approach,  the  Arrhenius  pre-factor  of  the  cesium  forward 
reaction  mechanisms,  A,  was  increased  by  factors  of  10  and  100,  respectively.  Con¬ 
ductivity  and  Schlieren-type  plots  are  shown  in  Figures  7.23  and  7.24  for  the  gener¬ 
ator  and  accelerator  configurations,  respectively.  One  can  see  a  significant  increase 
in  conductivity  of  the  detonation  close  to  the  leading  shock  for  both  the  accelera¬ 
tor  and  generator  configurations,  with  an  increase  in  the  “enhanced  kinetic  factor” , 
EK ,  where  the  modified  Arrhenius  pre-factor  in  Equation  2.24  is  A'  =  EK  x  A. 
Slight  deceleration  and  acceleration  of  the  detonation  front  is  observed,  respectively, 
in  Figures  7.23  &  7.24.  The  centerline  peak  pressure  traces  for  the  accelerator  config¬ 
uration,  illustrated  in  Figure  7.25,  show  that  as  the  EK  is  increased,  the  amplitude 
of  the  peak  pressure  decreases.  In  contrast,  the  centerline  peak  pressure  traces  for 
the  generator  configuration,  illustrated  in  Figure  7.26,  show  that  as  EK  is  increased 
the  amplitude  of  the  peak  pressure  increases.  These  results  are  consistent  with  the 
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ID  results  from  Section  7.2  for  the  accelerator  and  generator  configurations,  in  that 
when  there  is  sufficient  conductivity,  with  a  larger  concentration  of  Cs  for  the  ID 
detonation,  the  applied  magnetic  and  electric  fields  affect  the  amplitude  of  the  peak 
pressure.  But  even  with  the  large  increase  in  conductivity  shown  in  Figures  7.23  and 
7.24  for  EK  =  100,  the  X-t  diagram  of  the  unmagnetized,  accelerator,  and  gener¬ 
ator  configurations,  illustrated  in  Figure  7.27,  shows  only  a  small  alteration  of  the 
detonation  velocity:  +30m/s  (increase)  and  —20 m/s  (decrease)  for  the  accelerator 
and  generator  configurations,  respectively. 

7.3.2  Strong  Applied  Fields 

Next,  we  subjected  the  2D  detonation  to  significantly  stronger  magnetic  fields  (16T 
and  32  T)  without  enhancing  the  cesium  kinetics.1  Figures  7.28  and  7.29  show 
the  centerline  peak  pressure  traces  for  Bz  =  8,  16,  and  32  T  in  the  accelerator  and 
generator  configurations,  respectively.  Even  when  the  2D  detonation  in  the  generator 
configuration  is  subjected  to  a  strong  field,  the  peak  pressure  amplitude,  shown  in 
Figure  7.29,  is  slightly  increased.  The  X-t  plot  for  the  generator  configuration, 
illustrated  in  Figure  7.30,  shows  only  a  marginal  decrease  in  the  detonation  velocity. 
By  comparing  the  profiles  of  the  x-velocity,  conductivity,  pressure,  and  temperature 
for  the  generator,  illustrated  in  Figures  7.31  -  7.34,  respectively,  one  will  see  that  the 
MHD  effects  are  manifested  significantly  far  downstream  of  (>>  Aw)  the  leading 
shock. 

The  accelerator  configuration  with  a  strong  field,  on  the  other  hand,  has  a  more 

significant  effect  on  the  2D  detonation  dynamics.  From  Figure  7.28  one  can  see 

a  transition,  from  Bz  =  8 T  where  the  amplitude  and  frequency  of  oscillations  are 
1  Unrealistic  but  only  interested  in  scaling  and  dynamics 
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regular,  to  Bz  =  16 T  where  the  amplitude  begins  to  noticeably  decrease  at  f  ~  50/is, 
and  finally  to  Bz  =  32 T  where  at  t  =  40 /is  there  is  a  significant  decrease  in  amplitude 
and  an  increase  in  the  mean  peak  pressure.  Because  of  the  significant  change  in 
the  oscillation  frequency,  amplitude  and  mean  of  the  centerline  peak  pressure,  the 
associated  smoke-foil  record,  illustrated  in  Figures  7.35  -  7.38,  for  field  strengths 
Bz  =  0,  8T,  16T  and  32T,  respectively,  can  be  used  to  see  the  alteration  of  the 
underlying  structure  of  the  cellular  detonation.  In  Figure  7.36,  we  see  that  for  the 
Bz  =  8 T  case,  cellular  structure  pattern  remains  regular,  and  from  Figure  7.37  for  the 
Bz  =  16T  case,  the  cellular  patterns  become  irregular  at  x  ~  9.5 cm.  The  detonation 
case  where  Bz  =  32T  (Figure  7.38)  shows  that  there  is  a  transition  from  a  cellular 
structure  to  quasi-ID  detonation  at  x  ~  7.2 cm.  The  X-t  plot  for  the  accelerator 
configuration,  illustrated  in  Figure  7.39,  shows  marginal  increases  in  the  detonation 
velocity  of  +5 m/s  for  the  B ,  =  16,  but  there  is  an  increase  of  338 m/s  for  Bz  =  32 T, 
a  21%  increase  above  the  non-MHD  case. 

7.4  Conclusions 

In  this  chapter,  we  examined  the  effects  of  a  diluent  (cesium)  on  the  dynamics  of  a 
ID  spark-ignited  detonation  and  confirmed  the  observation  of  Radulescu  et  ah  [19] 
that  the  diluent  had  a  regularizing  effect  on  the  oscillations  of  the  ID  detonation. 
We  also  studied  a  cesium-seeded  ID  spark-ignited  detonation  subjected  to  an  ap¬ 
plied  magnetic  field  in  both  a  accelerator  and  generator  configurations  and  found 
that  the  accelerator  mode  had  a  regularizing  effect  on  the  detonation  oscillations, 
while  the  generator  mode  had  the  opposite  effect  (see  Figure  7.9).  While  it  was 
demonstrated  that  the  dynamics  of  a  cesium-seeded  ID  spark-ignited  detonation  can 
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be  significantly  altered  when  subjected  to  an  applied  magnetic  field  of  reasonable 
strength  (Bz  <  8 T),  the  same  statement  does  not  hold  true  for  the  two-dimensional 
case.  When  the  enhanced  kinetics  factor  EK  for  a  cesium-seeded  2D  spark-ignited 
detonation  was  increased  (EK  =  1,10,100),  the  conductivity  behind  the  leading 
shock  was  significantly  increased  (see  Figures  7.23  and  7.24),  but  there  was  no  alter¬ 
ation  to  the  detonation  velocity  (see  Figures  7.30  and  7.31).  The  2D  spark-ignited 
detonation  required  a  significantly  stronger  magnetic  field  (32 T)  than  the  ID  deto¬ 
nation  to  accelerate  the  detonation  in  the  time  (75 /is)  and  length  (15cm)  scales  of 
the  problem  (see  Figure  7.39).  These  findings  suggest  that  multidimensional  effects 
play  an  important  role  in  MHD  acceleration.  Perhaps  a  set  of  reactants/diluent  with 
an  increased  flame  temperature,  for  increased  conductivity,  and  a  lower  density,  for 
less  inertia,  would  be  more  amenable  for  MHD  acceleration. 
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%Cs 

Theoretical 
CJ  [m/s] 

1-D  [m/s] 

2-D  [m/s] 

0 

1967 

1998 

1938 

1 

1917 

1950 

1883 

5 

1754 

1763 

1714 

10 

1584 

1600 

1551 

Table  7.1:  The  effects  of  the  addition  of  Cesium  on  the  detonation  velocity  of  a 
ID  and  2D  H2  —  Air  detonation.  The  theoretical  Chapman- Jouguet(CJ)  velocity  is 
calculated  with  the  initial  conditions:  Pq  =  1  atm  and  T0  =  300A7 
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Figure  7.1:  Peak  pressure  traces  of  spark- ignited  H2-a ir  detonations  with  different 
amounts  of  seeded  cesium  without  an  applied  magnetic  field  (B  =  0).  The  MP5 
scheme  with  Ax  =  5 /zra  was  used  here. 


159 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


Figure  7.2:  Configuration  of  spark-ignited  detonation  with  an  applied  magnetic  field. 
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Figure  7.3:  Peak  pressure  traces  of  detonations  seeded  with  1%  cesium  subjected  to 
various  magnetic  field  strengths  Bz  without  an  applied  electric  field  (K  =  0).  The 
MP5  scheme  with  Ax  =  5 /im  was  used  here. 


161 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


(a)  Bz  =  OT 


Figure  7.4:  Peak  pressure  traces  of  detonations  seeded  with  5%  cesium  subjected  to 
various  magnetic  field  strengths,  Bz  without  an  applied  electric  field  (K  =  0).  The 
MP5  scheme  with  Ax  =  5 fj,m  was  used  here. 


162 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


60 


50 


- 

Oh 


20 


20  40  60  80  100  120  140 

time  [/is] 

(c)  Bz  =  5 T 


(d)  Bz  =  6 T 


Figure  7.5:  Peak  pressure  traces  of  detonations  seeded  with  10%  cesium  subjected 
to  various  magnetic  field  strengths  Bz  without  an  applied  electric  field  (K  =  0).  The 
MP5  scheme  with  Ax  =  5 /im  was  used  here. 


163 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


Figure  7.6:  Peak  pressure  traces  of  detonations  seeded  with  1%  cesium  subjected  to 
various  magnetic  field  strengths  Bz  with  an  applied  electric  field  (Ky  =  1.5),  for  an 
“accelerator”  configuration.  The  MP5  scheme  with  Ax  =  5 fim  was  used  here. 
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Figure  7.7:  Peak  pressure  traces  of  detonations  seeded  with  5%  cesium  subjected  to 
various  magnetic  field  strengths  Bz  with  an  applied  electric  field  (Ky  =  1.5),  for  an 
“accelerator”  configuration.  The  MP5  scheme  with  Ax  =  5 fim  was  used  here. 
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Figure  7.8:  Peak  pressure  traces  of  detonations  seeded  with  10%  cesium  subjected 
to  various  magnetic  field  strengths  Bz  with  an  applied  electric  field  ( Ky  =  1.5),  for 
an  “accelerator”  configuration.  The  MP5  scheme  with  Ax  =  5 fim  was  used  here. 
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(a)  Accelerator  Peak  pressure  trace  (b)  Generator  Peak  pressure  trace 


(c)  Accelerator  Induction  Length  (d)  Generator  Induction  Length 


(e)  Accelerator  Detonation  Velocity  (f)  Generator  Detonation  Velocity 


Figure  7.9:  Time  trace  comparisons  of  Bz  =  0,  3,  and  8 T  ID  detonations  with  10% 


Cs  in  the  accelerator  configuration {Ky  =  1.5)  and  generator  configuration (Ky  =  0.5) 


with  Ax  =  5/im. 
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Figure  7.10:  Computational  setup  for  2D  detonation  simulations,  where  the  red 
region  represents  the  computational  spark. 
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(a)  t  =  71.16 /rs 


(c)  t  =  71.64 ns 


(e)  t  =  72.12/xs 


(b)  t  =  71.4/zs 


(d)  t  =  71.88ns 


(f)  t.  =  72.36/is 


Figure  7.11:  Schlieren-type  plot  using  density  gradients  of  the  detonation  front  from 
time  71.16 fxs  to  72.36 /is. 
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(b)  Maximum  pressure  contour  with  linear  color  map  ( min  max) 


Figure  7.12:  Stoichiometric  H2—Air  detonation  at  t  =  75 /is  where  Ax  =  Ay  =  50 /xm 
and  x  :  y  e  [12, 15] cm  :  [0, 0.5]cm. 
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(a)  Conductivity  distribution  with  a  linear  color  map  overlaid  on  Schlieren-type  plot  using  density 
gradients  of  the  detonation  front 
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(b)  Maximum  pressure  contour  with  linear  color  map  ( min  max) 


Figure  7.13:  Stoichiometric  H2  —  Air  —  1  %Cs  detonation  at  t  —  75 /is  where 
Ax  =  Ay  =  50 ym  and  x  :  y  G  [11.5, 14.5]cm  :  [0,  0.5]cm. 
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(a)  Conductivity  distribution  with  a  linear  color  map  overlaid  on  Schlieren-type  plot  using  density 
gradients  of  the  detonation  front 
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(b)  Maximum  pressure  contour  with  linear  color  map  ( min  max) 


Figure  7.14:  Stoichiometric  H2  —  Air  —  5 %Cs  detonation  at  t  =  75 /is  where 
Ax  =  Ay  =  50 ym  and  x  :  y  G  [10, 13]cm  :  [0,  0.5]cm. 
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(a)  Conductivity  distribution  with  a  linear  color  map  overlaid  on  Schlieren-type  plot  using  density 
gradients  of  the  detonation  front 
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(b)  Maximum  pressure  contour  with  linear  color  map  ( min  max ) 

Figure  7.15:  Stoichiometric  H2  —  Air  —  10%  Cs  detonation  at  t  =  75/rs  where 
Ax  =  Ay  =  50 ym  and  x  :  y  G  [9, 12]cm  :  [0,  0.5]cm. 
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(a)  Wide  view 


(b)  Enlarged  view 


Figure  7.16:  X-t  plot  comparing  the  progression  of  the  leading  shocks  at  the  centerline 
of  the  2D  detonation  with  no  MHD,  generator,  and  accelerator  configurations.  Here 
the  mixture  has  1%  Cs  and  Bz  =  ST. 
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(a)  Wide  view 


(b)  Enlarged  view 


Figure  7.17:  X-t  plot  comparing  the  progression  of  the  leading  shocks  at  the  centerline 
of  the  2D  detonation  with  no  MHD,  generator,  and  accelerator  configurations.  Here 
the  mixture  has  5%  Cs  and  Bz  =  ST. 
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80 


(a)  Wide  view 


(b)  Enlarged  view 

Figure  7.18:  X-t  plot  comparing  the  progression  of  the  leading  shocks  at  the  centerline 
of  the  2D  detonation  with  no  MHD,  generator,  and  accelerator  configurations.  Here 
the  mixture  has  10%  Cs  and  Bz  =  8 T. 
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(a)  Wide  view 


(b)  Enlarged  view 

Figure  7.19:  Peak  pressure  trace  of  the  centerline  of  the  leading  shock  of  the  2D  det¬ 
onation  with  no  MHD,  generator,  and  accelerator  configurations.  Here  the  mixture 
has  1%  Cs  and  Bz  =  8  T. 
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(a)  Wide  view 


(b)  Enlarged  view 

Figure  7.20:  Peak  pressure  trace  of  the  centerline  of  the  leading  shock  of  the  2D  det¬ 
onation  with  no  MHD,  generator,  and  accelerator  configurations.  Here  the  mixture 
has  5%  Cs  and  Bz  =  8 T. 
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Figure  7.21:  Peak  pressure  trace  of  the  centerline  of  the  leading  shock  of  the  2D  det¬ 
onation  with  no  MHD,  generator,  and  accelerator  configurations.  Here  the  mixture 
has  10%  Cs  and  Bz  =  8 T. 
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Figure  7.22:  Centerline  pressure  and  conductivity  profiles  for  the  accelerator  and 
generator  configurations  with  Bz  =  8 T  at  t  =  75 [is. 
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Figure  7.23:  Conductivity  distribution  with  a  linear  color  map  overlaid  with 
Schlieren-type  plot  using  density  gradients  of  the  detonation  front  in  the  genera¬ 
tor  configuration  with  10%  Cs  and  Bz  =  8 T,  for  different  enhanced  kinetic(i7A") 
factors. 
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Figure  7.24:  Conductivity  distribution  with  a  linear  color  map  overlaid  with 
Schlieren-type  plot  using  density  gradients  of  the  detonation  front  in  the  acceler¬ 
ator  configuration  with  10%  Cs  and  Bz  =  8 T,  for  different  enhanced  kinetic(i7A") 
factors. 
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(a)  Wide  view 


(b)  Enlarged  view 


Figure  7.25:  Centerline  peak  pressure  trace  of  2D  detonation  for  the  accelerator 
configurations  with  10%  Cs  and  Bz  =  8 T  with  various  enhanced  kinetics  values  EK. 
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(a)  Wide  view 


(b)  Enlarged  view 


Figure  7.26:  Centerline  peak  pressure  trace  of  2D  detonation  for  the  generator  con¬ 
figurations  with  10%  Cs  and  Bz  =  8 T  with  various  enhanced  kinetics  values  E K . 
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(a)  Wide  view 


(b)  Enlarged  view 

Figure  7.27:  X-t  plot  of  2D  detonation  comparing  the  progression  of  the  leading 
shocks  at  the  centerline  of  the  no  MHD,  generator,  and  accelerator  configurations 
with  EK  =  100,  10%  Cs  and  Bz  =  8 T. 
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(a)  Wide  view 
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(b)  Enlarged  view 


Figure  7.28:  Centerline  peak  pressure  trace  of  2D  detonation  for  the  accelerator 


configuration  with  10%  Cs  and  various  Bz  values  without  enhanced  kinetics. 
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Figure  7.29:  Centerline  peak  pressure  trace  of  2D  detonation  for  the  generator  con¬ 
figuration  with  10%  Cs  and  various  Bz  values  without  enhanced  kinetics. 
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(a)  Wide  view 


(b)  Enlarged  view 


Figure  7.30:  X-t  plot  of  2D  detonation  comparing  the  progression  of  the  leading 
shocks  at  the  centerline  of  the  generator  configuration  with  10%  Cs  and  various  Bz 
values  without  enhanced  kinetics. 
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(a)  t  =  25 ns  (b)  t  =  37.5 fis 


(c)  t  =  50 /j,s  (d)  t  =  62.5 fis 


(e)  t  =  75 /j,s 


Figure  7.31:  X-velocity  profiles  of  2D  detonation  at  different  times  for  the  generator 


configuration  with  10%  Cs  and  various  Bz  values  without  enhanced  kinetics. 
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Figure  7.32:  Conductivity  profiles  of  2D  detonation  at  different  times  for  the  gener¬ 


ator  configuration  with  10%  Cs  and  various  Bz  values  without  enhanced  kinetics. 
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Figure  7.33:  Pressure  profiles  of  2D  detonation  at  different  times  for  the  generator 


configuration  with  10%  Cs  and  various  Bz  values  without  enhanced  kinetics. 
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(a)  t  =  25 ns  (b)  t  =  37.5 fis 


(c)  t  =  50 /j,s 


(d)  t  =  62.5 fis 


(e)  t  =  75 


Figure  7.34:  Temperature  profiles  of  2D  detonation  at  different  times  for  the  gener¬ 


ator  configuration  with  10%  Cs  and  various  Bz  values  without  enhanced  kinetics. 


192 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


-0.25 


0.0 


0.5 


-0.25 


0.0 


6.0 


0.0 


3.0 


P m ax  >  atm 
80.0 


70.0 


9.0 


60.0 


50.0 


40.0 


30.0 


20.0 


10.0 


9.0 


*  0.0 


12.0 


Figure  7.35:  Detonation  history  presented  by  maximum  pressure  contours  without 
MHD(5Z  =  OT)  for  10%  Cs. 
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Figure  7.36:  Detonation  history  presented  by  maximum  pressure  contours.  Acceler¬ 
ator  configuration  with  Bz  =  8 T  and  10%  Cs. 
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Figure  7.37:  Detonation  history  presented  by  maximum  pressure  contours.  Acceler¬ 
ator  configuration  with  Bz  =  16T  and  10%  Cs. 
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Figure  7.38:  Detonation  history  presented  by  maximum  pressure  contours.  Acceler¬ 
ator  configuration  with  Bz  =  32 T  and  10%  Cs. 
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(a)  Wide  view 


(b)  Enlarged  view 


Figure  7.39:  X-t  plot  of  2D  detonation  comparing  the  progression  of  the  leading 
shocks  at  the  centerline  of  the  accelerator  configuration  with  10%  Cs  and  various  Bz 
values  without  enhanced  kinetics. 
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CHAPTER  8 


Conclusions  and  Future  Work 

The  present  computational  studies  have  enabled  both  a  global  exploration  of  PDRIME 
concepts  for  propulsive  devices,  as  well  as  a  detailed  exploration  of  the  underlying 
physics  of  detonation-MHD  interactions. 

For  the  standard  PDRIME  utilizing  a  constant  magnetic  field,  there  was  marginal 
gains  in  performance  over  the  baseline  PDRE  configuration  using  simplified  mod¬ 
elling  approaches.  We  investigated  closed  and  open  loop  control  of  the  magnetic 
field  utilizing,  temporal  and  temperature  controllers,  respectively.  We  found  that 
these  methods  of  control  had  marginal  performance  gains,  at  best.  When  the  flight 
Mach  number  was  significantly  lowered  (M  <  5),  we  found  significant  performance 
increases  without  the  use  of  controllers. 

The  stability  and  dynamics  of  a  ID  spark-ignited  detonation  with  complex  kinet¬ 
ics  were  also  investigated.  We  developed  a  model  to  understand  and  interpret  the 
coupling  between  the  reactive  and  fluid  mechanical/acoustic  phenomena.  We  found 
that  the  frequency  in  oscillations  of  the  peak  pressure  trace  is  inversely  related  to 
the  time  it  takes  for  an  acoustic  wave  to  propagate  from  the  flame  to  the  leading 
shock  and  the  entropy  wave  generated  from  the  perturbed  shock  to  travel  back  to 
the  flame.  We  verified  the  model  by  finding  out  it  was  in  agreement  with  the  peak 
pressure  cycles  extracted  from  the  ID  complex  kinetics  detonation  simulation. 
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After  gaining  insight  into  the  stability  of  the  ID  detonation,  we  investigated  the 
stability  of  a  cesium-seeded  detonation.  We  confirmed  previous  findings  by  Radulescu 
et  al. [19]  that  the  diluent  had  a  regularizing  effect  on  the  ID  detonation.  We  then 
applied  MHD  using  a  loading  factor,  Ky,  to  the  ID  cesium-seeded  detonation  with 
both  the  generator  and  accelerator  configurations  at  various  magnetic  field  strengths, 
BZl  and  concentrations  of  cesium.  We  found  that  for  a  given  concentration  of  cesium 
and  applied  field  strength  the  accelerator  had  a  regularizing  effect,  while  the  genera¬ 
tor  had  the  opposite  effect.  We  then  investigated  the  effect  the  applied  fields  had  on 
the  2D  cesium-seeded  spark-ignited  detonation  with  the  same  parameters  used  for 
the  ID  simulations  (i.e.,  Bz  and  cesium  concentration)  and  found  that  the  MHD  had 
little  to  no  effect  on  the  detonation  dynamics,  i.e.,  cellular  structure  and  detonation 
velocity.  Next,  we  increased  the  Arrhenius  pre-factor  of  the  cesium  forward  reaction 
mechanism  by  factors  of  10  and  100.  We  observed  a  significant  increase  of  conductiv¬ 
ity  near  the  leading  shock  of  the  2D  detonation,  but  little  change  in  the  detonation 
velocity  in  either  the  accelerator  or  generator  configurations.  Lastly,  we  significantly 
increased  the  strength  of  the  magnetic  field  (16T  and  32T)  without  enhancing  the 
cesium  kinetics.  The  generator  mode  had  little  effect  on  the  detonation  dynamics, 
while  the  accelerator  mode  significantly  altered  the  detonation  velocity. 

In  Chapter  7,  the  conductivity  of  the  fluid  played  an  important  role  in  the  MHD 
acceleration  of  the  detonation.  But  in  order  to  properly  calculate  the  conductivity, 
the  number  density  of  electrons  as  well  as  the  electron  temperature,  Te,  must  be 
determined  (per  Equation  2.27).  Section  2.3.3  describes  how  the  electron  energy  is 
coupled  to  the  bulk  fluid  via  a  two-temperature  model.  One  future  direction  of  this 
research  is  to  investigate  the  dynamics  of  the  detonation  in  ID  and  2D  simulations 
using  the  two-temperature  model  to  more  accurately  characterize  the  thermal  non- 
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equilibrium  effects  of  the  bulk  fluid.  Additionally,  the  kinetics  of  excitation  and 
ionization  of  the  gas  can  be  further  investigated  by  means  of  a  collisional  radiative 
model[75,  76],  commonly  used  in  the  study  of  gas  discharges [82,  83],  which  allow  for 
a  more  accurate  representation  of  the  atomic  state  distribution  function. 

The  simulations  performed  in  this  dissertation  utilized  high-order  accurate  meth¬ 
ods  which  significantly  reduced  numerical  dissipation.  The  high-order  accuracy  al¬ 
lowed  for  the  capturing  of  sharp  flow  features  (i.e.,  shocks  and  contact  disconti¬ 
nuities).  Oran  et  al.[84]  demonstrated  that  the  viscous  effects  did  not  change  the 
dynamics  of  the  leading  shock,  thus  did  not  have  a  significant  effect  on  the  detona¬ 
tion  cellular  structure.  But  unlike,  previous  works[68,  84,  85]  with  detonations,  the 
downstream  effects  play  an  important  role  in  the  overall  dynamics  of  the  detonation 
with  MHD  acceleration.  Future  studies  involving  full  ionization  kinetics,  including 
collisional-radiative  processes,  will  be  used  to  examine  these  processes  in  further 
detail.  In  addition,  a  more  physically  accurate  model  based  on  the  Navier-Stokes 
equations  and  including  species  diffusion  is  needed  to  correctly  model  the  transport 
of  the  fluid.  In  addition,  precursor  effects [86]  might  play  a  role  in  heating  the  up¬ 
stream  fluid,  which  then  alters  the  flame  temperature.  It  might  also  be  interesting  to 
investigate  the  thermal  losses  to  the  wall  as  well  as  the  interactions  of  the  boundary 
layer  with  the  shock[87]. 
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APPENDIX  A 


Reaction  Mechanism 


A.l  #2- Air  Reaction  Mechanism 

The  772-air  reaction  mechanism  and  Arrenhius  coefficients  used  in  the  present  study 
contains  19  reversible  elementary  mechanisms  composed  of  9  species[37]. 


Elementary  Mech. 

A 

T] 

Ea 

1-) 

H  +  02  ^ 0  +  077 

2.60E14 

0.000 

8400 

2.) 

0  +  H2  ^ 77  +  077 

1.80E10 

1.000 

4450 

3.) 

H2  +  OH  ^h  +  h2o 

2.20E13 

0.000 

2575 

4.) 

20 77  ^  H20  +  0 

6.30E12 

0.000 

545 

5.) 

H  +  OH  +  M  ^  H20  +  M 

2.20E22 

-2.00 

0.000 

6.) 

2  77 +  M  ^  772  +  M 

6.40E17 

-1.000 

0.000 

7.) 

H+O+M  ^ OH+M 

6.00E16 

-0.600 

0.000 

8.) 

20  +  M  —  02  +  M 

6.00E13 

0.000 

-900 

9.) 

H2  +  O2  ^ HO2  +  H 

1.00E14 

0.000 

28000 

10.) 

H  +  O2  +  M  ^ HO2  +  M 

2.107715 

0.000 

-500 

11.) 

H  +  HO2  ^2  OH 

1.407714 

0.000 

540 

12.) 

H  +  HO2  ^ 0  +  H20 

1.007713 

0.000 

540 

13.) 

O  +  HO2  ^  02  +  OH 

1.507713 

0.000 

475 
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14.) 

oh  +  ho2  ^±o2  +  h2o 

8.007712 

0.000 

0.000 

15.) 

H202  +  M  ^2 OH  +  M 

1.207717 

0.000 

22750 

16.) 

2H02  +  M  ^  H202  +  02 

2.007712 

0.000 

0.000 

17.) 

h  +  h2o2  ^  h2  +  ho2 

1.407712 

0.000 

1800 

18.) 

0  +  H202  ^  OH  +  H02 

1.407713 

0.000 

3200 

19.) 

oh  +  h2o2  ^  h2o  +  ho2 

6.107712 

0.000 

715 

A. 2  Cesium  Reaction  Mechanism 


The  cesium  reaction  mechanism  and  Arrenhius  coefficients  used  in  the  present  study 
contains  2  reversible  elementary  mechanisms  composed  of  Cs,  Cs+,  and  e~. 


Elementary  Mech. 

A 

V 

Ea 

!•) 

Cs  +  e~  Cs+  +  e~  +  e~ 

2.487714 

0.500 

45900 

2.) 

Cs  +  M  ^  Cs  +  e“  +  M 

2.487711 

0.500 

45900 
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APPENDIX  B 


Eigenvector  Matrices 


The  following  was  derived  for  a  multi-species,  two-temperature  eigen-system [35,  88]. 

B.l  Governing  Equation 

The  governing  equation  for  a  non-reactive  multi-species,  multi-temperature  three- 
dimensional  flow  without  LHS  is: 


Q  t  +  VnF  —  0 


(B.l) 


Here  the  vector  containing  the  conserved  variables,  Q,  and  the  flux  normal  to  control 
volume  surface  ,  F.  are: 


Q  = 


/ps\ 
pu 
B 

E* 

V  ^  J 


( 


F  = 


Ps^n 

puun  +  P*  n  -  ±B„  B 


\ 


Un  B  -  U Bn 

(E*  +  P*)un  -  ±Bn u  •  B 


(B.2) 
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B.2  Roe  Averaged  Weighting 


Before  the  eigensystem  can  be  determined,  the  Roe  Averaged  variables [46],  i.e. ,  mass 
fraction,  spatial  components  of  velocity  and  enthalphy,  must  first  be  determined  on 
a  flux  interface,  i  +  \.  First  the  primitive  variables  are  calculated  for  i  and  i  +  1  as: 


/ 


Cs 


\ 


V  = 


u 

b 

h 


(B.3) 


\§e  J 

I  /  R2  \  n  9 

h  =  -  [E  +  P-- - Ee  )  ,  cs  =  — ,  b  =  [bn,  bt]T /yfpjT0  ,  and  se  =  — 

P  \  2p0  /  P  P 

where  cs  represents  the  mass  fraction  of  the  sth  species,  h  represents  the  specific 

enthalpy,  and  se  represents  the  specific  electron  entropy.  Next,  by  using  the  densities 

at  i  and  i  +  1,  the  Roe-averaged  variables  at  the  face,  i  +  |,  are  determined  by: 


i\fPi  4“  yj Pi-lf- 1 

yj Pi  yj  Pi+ 1 


(B.4) 


B.3  Eigensystem  and  Flux  Jacobian  Matrix 

The  flux  Jacobian  represents  the  relation  A  —  From  here  we  diagonalize  the 
flux  Jacobian  to  get  a  formulation  involving  the  Jacobian,  left  and  right  eigenvector 
matrices,  and  the  eigen-matrix. 

i=fh=RAL  (B.5) 
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The  eigen-matrix  is 

defined 

as  the  following, 

^  Un 

0 

0 

0 

0 

0 

0 

0 

0 

0  ^ 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0  xin 

+  cf 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0  u 

n  —  cf 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

U n  ~b  Cs 

0 

0 

0 

0 

0 

A  = 

0 

0 

0 

0 

0 

Un 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

un  cs 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Un  T  Ca 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Un 

0 

l 0 

0 

0 

0 

0 

0 

0 

0 

0 

Un  / 

(B.6) 

Where  the  diagonal  entries  are  the  eigen-values  of  the  system.  The  definitions 
of  the  fast(slow)  magneto-acoustic,  c/(s),  and  the  Alfven,  ca,  wave  speeds  can  be 
found  in  [88].  The  similarity  transformation  matrices  R  and  L  are  defined  in  the 
next  sections. 


B.3.1  Right  Eigenvectors 

The  system  satisfies  A Q  =  R  •  a  where 
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R  = 


1 

rf 

rf 

rs 

0 

rs 
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0 

'Cj x 

®t,x 
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y/~Psx 

y/~Psx 

uy 

Qtv 

®f,y 

et„ 
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yfpSy 

y[Psy 

uz 
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n7 
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0 

0 

0 

0 

0 

0 

ue 

1  ) 


(B.7) 


where  sB  is  the  sign  of  the  normal  component  of  the  magnetic  held  and  for  the 
purpose  of  a  compact  expression  of  R,  the  following  variables  are  defined, 


(B.8a) 


~  a 2  Cf 

A.  =  -r/-2 

c  f  yfp 

f  /  =  ~sBbnrs 


r  s  =  +sBarf 


ue  = 


Ce 


7e 


-P 


i7e-i 


(B.8b) 

(B.8c) 

(B.8d) 

(B.8e) 


~  B  ~ 

=  rf(s)h  ±  rmUnCf{s)  ±  Ut f/W  +  ^A/(s)  (B.8f) 

w  ho 

0/(s),x  =  rf(s)Ux  A  r  f{s)cf{s)nx  i  f  (B.8g) 

similarly  for  the  y,  z  components  and  £e  =  1  —  As  noted  by  Brio  &  Wu[60], 
renormalization  factors,  r/  and  rs,  are  needed  to  avoid  singular  solutions  when  the 
magnetic  held  vanishes.  Their  dehnitions  and  identities  are  as  follows1, 

i  b2  r 

Identities  computed  with  et  =  a2+b2  >  ~  10“ 
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(B.9a) 


B.3.2  Left  Eigenvectors 

The  systems  satisfies  L  =  R 


(B.9b) 


(B.9c) 


2,  2^  _  1 

rs+rf  r2  -  1 
°f 


(B.9d) 


r/  +  rH  =  1 

J  az 


(B.9e) 


1 ,  where 


/ 


\ 


Lf 


Lf 

Lpseudo 

Li 

L\ 

L~a 

Le 
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—  SeP\ 
f3uxSe 
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B.3.3  Riemann  ‘Jump’  conditions 


The  ‘jump’  conditions  satisfy  a  =  L  ■  AQ  and  is  defined, 


/A 


a. 


a. 


a  = 


« 8 
a~ 

dj 

a 


A 


where, 


Oir 


=  1 


V  / 


ph) Ap + 5A(pS) + ^ABi  ~  !aae‘ 


(B.ll) 


(B.12a) 
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(B.12b) 


[AP  ±  pcfAun ]  +rs^~ 
Zcf 


A Bt  =f  sByfp—Aut 
cf 


[AP  ±  pcsAun]  -  rf^~  ABtT  ssy/p—Aut 

/Cf  L  CL 


(B.12c) 


a  A  = 


1 

2 


[y/pA us  =F  sbABs] 


(B.12d) 


0Le 


A  se  - 


EeAP 

p  a 2 


(B.12e) 
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APPENDIX  C 


MHD  Divergence  Cleaning  for  General 
Coordinate  Systems 


The  pseudo  electric  field  is  defines  as  hi  =  u  x  B.  From  Maxwell’s  Equations  the 
temporal  evolution  can  be  described  by 


<9B 

~dt 


-V  x  12 


(C.l) 


For  a  cartesian  coordinate  system  with  x,  y,  and  z  dependants  ^  is  described  as: 


9BX 

—  —12  - 

dy  * 

e_t 

dt 

dz' 

9By 

—  —12  - 

dz^x 

_d_ 

dt 

dx 

9BZ 

—  —12  - 

_d_ 

dt 

9y 

(C.2) 


For  a  coordinate  system  with  a  r  —  z  dependants  ^  is  described  as: 


°-§r  =-iii(rne) 

B~§?  =l&(rQo) 

9Bg  _ _ do  |  d_r\ 

dt  —  9rlizTaz 

For  a  coordinate  system  with  a  r  —  9  dependants  is  described  as: 


(C.3) 


9Br 

—  4-5-12 

r  d6»“2 

9t 

9BZ 

1  d  ( 

9t 

r  ae 

9Be 

=  -—12 
dr 

9t 

r  dr* 


(C.4) 
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APPENDIX  D 


Jacobians  and  Transforms 


D.l  Chemical  Jacobian 


We  first  start  out  with  the  model  ordinary  differential  equation  for  the  kinetics: 

where 

QT')  *-(£) 

ns  is  the  number  density  of  the  sth  species,  and  Cjs  and  uje  are  the  species  and 
energy  production  terms,  respectively.  The  number  density,  ns,  was  used  in  lieu 
of  mass  density,  ps,  out  of  convience  and  can  be  easily  transformed,  =  MkSki- 
The  implicit  1st  order  numerical  formulation  of  Equation  D.l  is  carried  out  with  the 
following  steps: 

AQ  _  An+1 
At  “ 

AQ  (S n  I  dSi  A  f 

At  +  9t“ 

AQ  _  nn  ,  t)il  AQ  a  + 

At  +  3Q  A t  ^ 

%  =sWf§AQ 


(D.3) 


I  -  IgA t )  AQ  =  ft"  A t 


AQ  =  ( I  —  §§At)  1  ftnAt 
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D.2  Chemical  Jacobian(|^)  Derivation 


^  =  TT n'rj 

dnk  ^rsnkH> 


dcoq  du.  dT  du.  1 


d E  dT  dE  dT  Cv 
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^  rs  dT  11  k 
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r  'Ir 


—ArTVrexp 


0r  \  0r 


T  +  —ArTVrexp  T 
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dE  Cv 
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rp  p  2 
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du  e  v  ■>  du  s 
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Quje  v  t  dio  s 
~dE  =  ^  ~dEe°s 

S 


(D.4) 


(D.5) 


(D.6) 


(D.7) 
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APPENDIX  E 


Steady  State  Detonation 

We  will  first  begin  with  our  model  equation  with  the  steady  state  approximation: 

(E.la) 

(E.lb) 
(E-lc) 


From  Equation  E.lc,  the  explicit  space-marching  formulation  is  defined  as: 

^9  A-'a 

Ax 

The  implicit  space-marching  formulation  is  also  defined, 


(E.2) 


1*  -  A^+‘ 

(E.3a) 

(E.3b) 

(E.3c) 

Now  by  solving  for  AQ  the  final  form  is  expressed  as: 

AQ  =  (  A  -  Ax^^j  AAx 


(E.4) 
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APPENDIX  F 


Iterative  &  Direct  Solvers 


F.l  Thomas’  Algorithm 


The  triadiagonal  matrix  algorithm,  commonly  refered  to  as  Thomas’  algorithm,  is  a 
simplified  form  of  Gaussian  elimination  used  to  solve  triadiagonal  system  of  equation. 
These  systems  of  equations  take  on  the  following  form: 


(F.l) 


where  a o  =  0  and  cn  =  0.  And  represented  in  the  matrix  form  as: 


1 

o 

<o 

-o 

X\ 

d\ 

02  C-2 

X2 

ft- 

to 

a3  b3  . 

= 

cn-  i 

0  Ojv  b n 

.  Xn  . 

.  d'  N  . 

(F.2) 


Just  as  in  Gaussian  elimination,  Thomas’  algorithm  consist  of  a  forward  elimination 
and  backward  substitution  to  solve  the  system  as  follows: 

Forward  Elimination 
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for  k  =  2  loop  through  N 


Ik 

bk- 1 

bk  TYlCk—1 

dk  -  mdk- 1 

dN 

*  ~  -j— 

bN 

ik  CkXk-\-l 
bk 

This  algorithm  is  applicable  for  diagonally  dominant  matrices,  where 

\bi\  >  \di\  +  \ci\  i  E  1, N  (F.3) 

F.2  Black-Red  Gauss-Seidel 

Gauss-Seidel  is  an  iterative  method  used  to  solve  problems  of  the  general  form: 

U  =  /( U)  (F.4) 

If  one  were  to  solve  this  equation  explicitly,  it  would  simply  take  the  form, 

U-+1  =  /( un)  (F.5) 

If  U  is  the  spatially  distribution  of  U,  Equation  F.4  can  be  solved  implicitly  using 
Gauss-Seidel.  In  this  procedure,  one  would  begin  with  an  initial  solution  at  t  =  tn+1, 
TjW,  where  s=0  initially.  The  solution  at  is  used  to  determine  U^+1^  until 
convergence  which  is  as  follows: 
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m  = 
bk  = 

dk 

Backward  Substitution 

xA 

for  k  =  N  —  1  loop  through  1 
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error  =  max 


(F.6) 


I  I  X  j(s+l)  T  j(s)  I  I 

continue  until  ... 

error  <  e  — »  Un+1  =  TTS+1)  (convergence!) 

Using  this  procedure,  the  update  solution,  U^+1^,  is  solved  for  using  the  latest  solu¬ 
tion.  In  order  to  speed  up  this  process,  instead  of  sweeping  through  all  of  the  cells 
(i  &  j  for  2D),  one  can  first  sweep  though  and  solve  for  the  computational  grid  cells 
colored  red  in  Figure  F.l,  then  using  the  updated  solutions,  sweep  through 

and  solve  for  the  black  grid  cells,  The  procedure  would  be  as  follows: 

t  j(s+1)  _  /Yfj(s)  ^ 
ure<2  —  J  {V  black) 

black  ~  J  \  red  ) 

error  =  max  I  |U-S.+1^  —  U,-s,- 1 1  (F.7) 

I  I  l,J  I  I 

continue  until  ... 

error  <  e  — »  Un+1  =  Lhs+1)  (convergence!) 

The  procedure  shown  above  is  refered  to  as  Red-Black  Gauss-Seidel. 
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Figure  F.l:  Computational  domain  split  into  red  and  black  computational  cells. 
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APPENDIX  G 


Message  Passing  Interface  (MPI)  Implementation 

G.l  Grid  Connectivity 

A  three  dimensional  computational  domain  can  be  broken  up  into  rectangular  cuboid 
subdomains;  each  containing  6  rectangular  faces.  A  face  can  connect  with  one  or  more 
faces  from  adjacent  domain(s),  which  is  illustrated  in  Figure  G.l.  This  illustration 
also  shows  that  domain  pO  is  connected  with  2  domains(pl  &  p2),  domain  pi  is 
connected  with  3  domains(pl,p2,p3),  and  domain  p3  is  connected  with  2  domains  (pi 
&  p2).  The  red  numbers  in  the  figure  represent  unique  domain-to-domain  connection 
of  which  there  are  5.  Domains  pO  &  p2  of  Figure  G.l  are  illustrated  with  there 
associated  ghost  layers  iin  Figure  G.2.  For  this  connection,  labeled  1  in  Figure  G.l, 
pO  would  pass  information  of  its  interior  cells  to  the  ghost  layer  of  pi  and  pi  would  in 
turn  information  of  its  interior  to  p2.  This  example  shows  that  for  each  connection, 
there  are  two  receive  buffers  required,  which  will  be  refered  to  as  memory  windows 
without  explanation  for  the  moment.  From  the  5  connections  in  Figure  G.l,  a  list  of 
10  memory  windows  and  its  associated  domains  are  created  in  Table  G.l. 
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n 


Figure  G.l:  The  computational  domain,  Q  is  decomposes  into  4  subdomains. 


pO 


Figure  G.2:  Domain-to- domain  connection  example  where  the  red  shaded  area  is  a 
pO  ghost  to  p2  physical  window  and  the  grey  shaded  area  is  a  p2  ghost  to  pO  physical 
window. 


window 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

local 

0 

2 

0 

1 

2 

1 

2 

3 

3 

1 

remote 

2 

0 

1 

0 

1 

2 

3 

2 

1 

3 

Table  G.l:  Generic  memory  window  list  generated  from  Figure  G.l  domain  connec¬ 
tions  where  ‘local’  is  the  domain  sending  its  physical  cell  data  to  the  ghost  cell  of 
the  ‘remote’  domain. 
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